Setup

Loading Packages

We will start by loading some useful packages. Using the library function, we can load these packages if they have previously been installed in your R. If they have not been previously installed, replace the library function with install.packages function and place the package name in quotation marks. That will install the package. Next, run the library function again and it will load the package. You only need to install a package once.

library(plyr) # Necessary for Rmisc
library(Rmisc) # Necessary for summarySE
library(dplyr) # Necessary for group_by
library(reshape) # Necessary for melt & cast
library(ggplot2) # Necessary for ggplots
library(grid) # Necessary for pushViewPort
library(ggpubr) # Necessary for ggarrange
library(conover.test) # Necessary for conover test
library(lmtest) # Necessary for homoscedasticity test.
library(lme4) # Generalized linear models and mixed effects.
library(lmerTest) # Necessary to get p-vals in glms
library(lubridate) # Necessary for monthiness? 
library(ggthemes) # Necessary for cool themes 
library(emmeans) # Display all pairwise comparisons
library(ggfortify) # Necessary for PCA. 
library(patchwork) # Necessary for mixing ggplots and legends
library(kableExtra) # Necessary for 'kbl' function.
library(report) # Necessary for 'sessionInfo' function.

Reading Dataframes

We will set the working directory using the setwd function. If you are looking to run this code in your computer, be sure to change the path to your data. load dataframes using the read.csv function.

setwd("~/Desktop/R") 
LeafLitt <- read.csv("Data/Leaf.Litter.Ecosystems.csv") 
Clim <- read.csv("Data/Clim.csv") 
Struc <- read.csv("Data/StrucDF.csv") 
Solar <- read.csv("Data/SolarRadiation.csv") 
Irad <- read.csv("Data/SolarIrradiation.csv") 
Composition <- read.csv("Data/PlotStructureCaqueta.csv") 

Now, the next step is to get these dataframes in the format that we want. Several columns are of the character class but we want them as factors. Here the across function selects the columns to be converted to factors based on their names, and as.factor is applied to each selected column using the mutate function from the dplyr package. Finally, the resulting factors are assigned back to the corresponding columns using the %>% pipe operator.

Correcting Variable Class

LeafLitt <- LeafLitt %>%
  dplyr::mutate(across(c("Paisaje", "Sample", "Locality", "Successional.State", "Succession",
                         "Plot", "Record"), as.factor))
Clim <- Clim %>%
  dplyr::mutate(across(c("Landscape", "Plot", "Sample"), as.factor))

Struc <- Struc %>%
  dplyr::mutate(across(c("Landscape", "Plot", "Succession", "CodeZ"), as.factor))

Solar <- Solar %>%
  dplyr::mutate(across(c("Estacion", "Municipio"), as.factor))

Irad$Month <- as.factor(Irad$Month) # Florencia Data. 

Composition <- Composition %>%
  dplyr::mutate(across(c("Landscape", "Chrono", "Chronosecuence", "Plot", "Species"), as.factor))

Great. Next we need to change some values from Muestro 5 (353 (DM2), 365(RM), and 377(RM)) into NAs, given that Muestreo 5 is incomplete for some plots.

Setting NAs

LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$LeafLitter.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Branches.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Detritus.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Flowers.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Fruits.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "353")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "365")] <- NA
LeafLitt$Total.Mg.ha.month[which(LeafLitt$Record == "377")] <- NA

Dataframe Setup

In the LeafLitt dataframe, there is one level of the variable successiont hat we will be removing, DV1. We will create a new dataframe that excludes this level using the droplevels and subset functions.

LeafLitt2 <- droplevels(subset(LeafLitt, LeafLitt$Succession != "DV1"))

Next, we will use the paste function to create a new variable that is a combination of the text from the Paisaje, Successional.State, and Plot (CodeX). We will also create one for a combination of Succession and Plot (CodeZ). These variables will be used for matching dataframes further down the code.

LeafLitt2$CodeX <- paste(LeafLitt2$Paisaje , LeafLitt2$Successional.State, LeafLitt2$Plot)
LeafLitt2$CodeZ <- paste(LeafLitt2$Succession , LeafLitt2$Plot)

Here we will create a CodeX for the Clim dataframe as well. Then we will use CodeX to attach the Succession column from LeafLitt2 using the with and match functions. Next we remove NA values from Clim dataframe using droplevels and subset functions to create a new dataframe, Clim2. A CodeY variable is created by pasting CodeX and Sample for Clim2 and LeafLitt2.

Clim$CodeX <- paste(Clim$Landscape , Clim$Successional.State, Clim$Plot)
Clim$Succession <- with(LeafLitt2, Succession[match(Clim$CodeX, CodeX)]) 
Clim2 <- droplevels(subset(Clim, Clim$Succession != "NA"))
Clim2$CodeY <- paste(Clim2$CodeX , Clim2$Sample)
LeafLitt2$CodeY <- paste(LeafLitt2$CodeX , LeafLitt2$Sample)

Using the shared CodeY, we will now insert the precipitation and temperature data into the LeafLitter2 dataframe. Note that Rain1 has a 1-mth lag, Rain2 has 2-mth lag, and Rain3 3-mth long. Same for temperature.

LeafLitt2$RainNow <- with(Clim2, RainNow[match(LeafLitt2$CodeY, Clim2$CodeY)]) 
LeafLitt2$Rain1 <- with(Clim2, Rain1[match(LeafLitt2$CodeY, Clim2$CodeY)]) 
LeafLitt2$Rain2 <- with(Clim2, Rain2[match(LeafLitt2$CodeY, Clim2$CodeY)]) 
LeafLitt2$Rain3 <- with(Clim2, Rain3[match(LeafLitt2$CodeY, Clim2$CodeY)]) 

LeafLitt2$TempNow <- with(Clim2, TempNow[match(LeafLitt2$CodeY, Clim2$CodeY)]) 
LeafLitt2$Temp1 <- with(Clim2, Temp1[match(LeafLitt2$CodeY, Clim2$CodeY)]) 
LeafLitt2$Temp2 <- with(Clim2, Temp2[match(LeafLitt2$CodeY, Clim2$CodeY)]) 
LeafLitt2$Temp3 <- with(Clim2, Temp3[match(LeafLitt2$CodeY, Clim2$CodeY)]) 

We’ll use the same trick as before to make our Code variables into factors.

LeafLitt2 <- LeafLitt2 %>%
  dplyr::mutate(across(c("CodeX", "CodeZ", "CodeY"), as.factor))

For various analyses we will want to deal with landscapes separately, so here we create one dataframe for Hill plots (Lomerio) and for Mountain using droplevels and subset functions.

Lomerio <- droplevels(subset(LeafLitt2, LeafLitt2$Paisaje == "Lomerio"))
Mountain <- droplevels(subset(LeafLitt2, LeafLitt2$Paisaje == "Mountain"))

Now, the months in these dataframes are coded as numbers, and these numbers differ by dataframe, so for each dataframe, we will manually relabel the months.

levels(Lomerio$Sample) <- list("Aug" = "1",  "Sept" = "2", "Oct"  = "3", 
                               "Nov" = "4", "Dec"  = "5", "Jan" = "6",
                               "Feb" = "7", "Mar" = "8", "Apr"  = "9", 
                               "May" = "10", "Jun"  = "11",  "Jul" = "12")

levels(Mountain$Sample) <- list("Sept" = "1",  "Oct" = "2", "Nov"  = "3", 
                                "Dec" = "4", "Jan"  = "5", "Feb" = "6",
                                "Mar" = "7", "Apr" = "8", "May"  = "9", 
                                "Jun" = "10", "Jul"  = "11",  "Aug" = "12")

We continue splitting, this time by Age Group.

LomeDL1 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL1"))
LomeDL2 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL2"))
LomeRL <- droplevels(subset(Lomerio, Lomerio$Succession == "RL"))

MounDM1 <- droplevels(subset(Mountain, Mountain$Succession == "DM1"))
MounDM2 <- droplevels(subset(Mountain, Mountain$Succession == "DM2"))
MounRL <- droplevels(subset(Mountain, Mountain$Succession == "RM"))

Descriptive Statistics

Mean and Standard Error

We can obtain mean and standard error for each variable in each landscape by running the summarySE function on the LeafLitt2 dataframe and setting Paisaje as the grouping variable. Before that, we will use the options function to set how many digits we want to see. We will wrap this up within the kbl function, which allows me to present the table in a neater format. We will also exclude some other columns using the square brackets.

options(digits=3)

kbl(summarySE(LeafLitt2, measurevar="LeafLitter.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Paisaje N LeafLitter.Mg.ha.month se
Lomerio 156 0.603 0.022
Mountain 213 0.433 0.015
kbl(summarySE(LeafLitt2, measurevar="Branches.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Paisaje N Branches.Mg.ha.month se
Lomerio 156 0.040 0.003
Mountain 213 0.051 0.003
kbl(summarySE(LeafLitt2, measurevar="Detritus.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Paisaje N Detritus.Mg.ha.month se
Lomerio 156 0.132 0.007
Mountain 213 0.135 0.006
kbl(summarySE(LeafLitt2, measurevar="Flowers.Mg.ha.month", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Paisaje N Flowers.Mg.ha.month se
Lomerio 156 0.002 0.000
Mountain 213 0.003 0.001
kbl(summarySE(LeafLitt2, measurevar="Fruits.Mg.ha.month", groupvars=c("Paisaje"), na.rm = 
TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Paisaje N Fruits.Mg.ha.month se
Lomerio 156 0.022 0.003
Mountain 213 0.018 0.002
kbl(summarySE(LeafLitt2, measurevar="Total.Mg.ha.month", groupvars=c("Paisaje"), na.rm = 
TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Paisaje N Total.Mg.ha.month se
Lomerio 156 0.807 0.028
Mountain 213 0.647 0.021

Percent Composition

Ok, so how much of the capture leaf litter is made of Leaf, Branch, Fruit, Flower, or Detritus?

options(digits=4)
DFX1 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Paisaje, data = LeafLitt2, FUN = sum, na.action=na.exclude)
DFX1$LeafP <- DFX1[,2]/DFX1[,7]*100
DFX1$BranchP <- DFX1[,3]/DFX1[,7]*100
DFX1$FlowerP <- DFX1[,4]/DFX1[,7]*100
DFX1$FruitP <- DFX1[,5]/DFX1[,7]*100
DFX1$DetritusP <- DFX1[,6]/DFX1[,7]*100
kbl(DFX1[c(1,8:12)], format = "html", table.attr = "style='width:50%;'") %>% kable_classic()
Paisaje LeafP BranchP FlowerP FruitP DetritusP
Lomerio 74.73 5.007 0.2344 2.717 16.38
Mountain 66.98 7.920 0.5166 2.800 20.93
DFX4 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Succession,
                  data = Lomerio, FUN = sum, na.action=na.exclude)

DFX4$LeafP <- DFX4[,2]/DFX4[,7]*100
DFX4$BranchP <- DFX4[,3]/DFX4[,7]*100
DFX4$FlowerP <- DFX4[,4]/DFX4[,7]*100
DFX4$FruitP <- DFX4[,5]/DFX4[,7]*100
DFX4$DetritusP <- DFX4[,6]/DFX4[,7]*100
kbl(DFX4[c(1,8:12)], format = "html", table.attr = "style='width:50%;'") %>% kable_classic()
Succession LeafP BranchP FlowerP FruitP DetritusP
DL1 77.01 3.654 0.0208 3.342 14.99
DL2 74.04 5.152 0.2038 3.072 16.61
RL 74.66 5.376 0.3598 2.017 16.66
DFX4 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Succession,
                  data = Mountain, FUN = sum, na.action=na.exclude)

DFX4$LeafP <- DFX4[,2]/DFX4[,7]*100
DFX4$BranchP <- DFX4[,3]/DFX4[,7]*100
DFX4$FlowerP <- DFX4[,4]/DFX4[,7]*100
DFX4$FruitP <- DFX4[,5]/DFX4[,7]*100
DFX4$DetritusP <- DFX4[,6]/DFX4[,7]*100
kbl(DFX4[c(1,8:12)], format = "html", table.attr = "style='width:50%;'") %>% kable_classic()
Succession LeafP BranchP FlowerP FruitP DetritusP
DM1 71.71 6.694 0.3250 1.651 18.70
DM2 64.04 8.865 0.4382 2.915 22.93
RM 66.74 7.797 0.7298 3.477 20.41

Diagnostic Statistics

Normality Test

We will use Shapiro-Wilk test, available through the shapiro.test function, to test for normality in our response variables. The Shapiro-Wilk test generates a test statistic, usually denoted as W. The closer the value of W is to 1, the more the data resemble a normal distribution. The shapiro.test function also provides a p-value, which represents the probability of obtaining the observed data if the null hypothesis (data are normally distributed) is true. If the p-value is greater than the chosen significance level (commonly 0.05), we fail to reject the null hypothesis, suggesting that the data can be considered approximately normally distributed. If the p-value is less than the chosen significance level, we reject the null hypothesis and conclude that the data significantly deviate from a normal distribution.

shapiro.test(LeafLitt2$LeafLitter.Mg.ha.month) 
## 
##  Shapiro-Wilk normality test
## 
## data:  LeafLitt2$LeafLitter.Mg.ha.month
## W = 0.94, p-value = 5e-11
shapiro.test(LeafLitt2$Branches.Mg.ha.month) 
## 
##  Shapiro-Wilk normality test
## 
## data:  LeafLitt2$Branches.Mg.ha.month
## W = 0.85, p-value <2e-16
shapiro.test(LeafLitt2$Detritus.Mg.ha.month) 
## 
##  Shapiro-Wilk normality test
## 
## data:  LeafLitt2$Detritus.Mg.ha.month
## W = 0.85, p-value <2e-16
shapiro.test(LeafLitt2$Flowers.Mg.ha.month) 
## 
##  Shapiro-Wilk normality test
## 
## data:  LeafLitt2$Flowers.Mg.ha.month
## W = 0.42, p-value <2e-16
shapiro.test(LeafLitt2$Fruits.Mg.ha.month) 
## 
##  Shapiro-Wilk normality test
## 
## data:  LeafLitt2$Fruits.Mg.ha.month
## W = 0.61, p-value <2e-16
shapiro.test(LeafLitt2$Total.Mg.ha.month) 
## 
##  Shapiro-Wilk normality test
## 
## data:  LeafLitt2$Total.Mg.ha.month
## W = 0.95, p-value = 2e-09

None of the variables follow a normal distribution. This isn’t necessarily wrong or bad, but it does mean that we may obtain more accurate models if we account for lack of normality.

Data Distributions

Let’s take a look at some distributions. For this we will create a histogram and kernel density estimate, which is basically a smoothed out histogram using the geom_histogram and geom_density arguments of ggplot. We will also use geom_vline to create a vertical line where the mean of our species richness is.

ggplot(LeafLitt2, aes(x=LeafLitter.Mg.ha.month)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
  geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
  geom_vline(aes(xintercept=mean(LeafLitter.Mg.ha.month, na.rm=TRUE)),color="cadetblue4", 
             linetype="dashed", size=1) + # This adds a blue line where the mean is. 
  xlab("Leaf Litter (Mg/Ha/Mth") + 
  ylab("Density") + 
  theme_bw(base_size = 15)

ggplot(LeafLitt2, aes(x=Branches.Mg.ha.month)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
  geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
  geom_vline(aes(xintercept=mean(Branches.Mg.ha.month, na.rm=TRUE)),color="cadetblue4", 
             linetype="dashed", size=1) + # This adds a blue line where the mean is. 
  xlab("Branch Litter (Mg/Ha/Mth") + 
  ylab("Density") + 
  theme_bw(base_size = 15)

ggplot(LeafLitt2, aes(x=Detritus.Mg.ha.month)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
  geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
  geom_vline(aes(xintercept=mean(Detritus.Mg.ha.month, na.rm=TRUE)),color="cadetblue4", 
             linetype="dashed", size=1) + # This adds a blue line where the mean is. 
  xlab("Detritus Litter (Mg/Ha/Mth") + 
  ylab("Density") + 
  theme_bw(base_size = 15)

ggplot(LeafLitt2, aes(x=Flowers.Mg.ha.month)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
  geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
  geom_vline(aes(xintercept=mean(Flowers.Mg.ha.month, na.rm=TRUE)),color="cadetblue4", 
             linetype="dashed", size=1) + # This adds a blue line where the mean is. 
  xlab("Flower Litter (Mg/Ha/Mth") + 
  ylab("Density") + 
  theme_bw(base_size = 15)

ggplot(LeafLitt2, aes(x=Fruits.Mg.ha.month)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
  geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
  geom_vline(aes(xintercept=mean(Fruits.Mg.ha.month, na.rm=TRUE)),color="cadetblue4", 
             linetype="dashed", size=1) + # This adds a blue line where the mean is. 
  xlab("Fruit Litter (Mg/Ha/Mth") + 
  ylab("Density") + 
  theme_bw(base_size = 15)

ggplot(LeafLitt2, aes(x=Total.Mg.ha.month)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="cornsilk") +
  geom_density(alpha=.2, fill="#FF6666") + # easier to see distribution
  geom_vline(aes(xintercept=mean(Total.Mg.ha.month, na.rm=TRUE)),color="cadetblue4", 
             linetype="dashed", size=1) + # This adds a blue line where the mean is. 
  xlab("Total Litter (Mg/Ha/Mth") + 
  ylab("Density") + 
  theme_bw(base_size = 15)

All distributions are right-skewed, some more than others.

Homoscedasticity & Linearity

Next let’s check for homoscedasticity and linearity. To do this, we’l create our first linear models.

lm1 <- lm(LeafLitter.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm1, 1) # Residual vs Fitted looks fine. 
plot(lm1, 3) # Scale-Location looks fine too. 

lm2 <- lm(LeafLitter.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm2, 1) # Residual vs Fitted looks fine. 
plot(lm2, 3) # Scale-Location looks fine too. 

lm3 <- lm(Branches.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm3, 1) # Residual vs Fitted looks fine. 
plot(lm3, 3) # Scale-Location shows increase of residuals with fitted values. 

lm4 <- lm(Branches.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm4, 1) # Residual vs Fitted looks fine. 
plot(lm4, 3) # Scale-Location looks fine too. 

lm5 <- lm(Flowers.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm5, 1) # Residual vs Fitted looks fine. 
plot(lm5, 3) # Scale-Location Shows mild increase.

lm6 <- lm(Flowers.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm6, 1) # Residual vs Fitted looks fine. 
plot(lm6, 3) # Scale-Location shows increase of residuals with fitted values. 

lm7 <- lm(Fruits.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm7, 1) # Residual vs Fitted looks fine. 
plot(lm7, 3) # Scale-Location looks fine too. 

lm8 <- lm(Fruits.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm8, 1) # Residual vs Fitted looks fine. 
plot(lm8, 3) # Scale-Location looks fine too. 

lm9 <- lm(Detritus.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm9, 1) # Residual vs Fitted looks fine. 
plot(lm9, 3) # Scale-Location looks fine too. 

lm10 <- lm(Detritus.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm10, 1) # Residual vs Fitted looks fine. 
plot(lm10, 3) # Scale-Location looks fine too. 

lm11 <- lm(Total.Mg.ha.month ~ Succession, data = Lomerio)
plot(lm11, 1) # Residual vs Fitted looks fine. 
plot(lm11, 3) # Scale-Location looks fine too. 

lm12 <- lm(Total.Mg.ha.month ~ Succession, data = Mountain,na.action=na.exclude)
plot(lm12, 1) # Residual vs Fitted looks fine. 
plot(lm12, 3) # Scale-Location looks fine too. 

Scedasticity and linearity are largely fine. We seem to mostly have a lack of non-normality with a right-skew.

Transformations

For left-skewed data—tail is on the left, negative skew, common transformations include (i) square root, (ii) cube root, and (iii) log.

First test: Square root transformation.

hist(sqrt(LeafLitt2$LeafLitter.Mg.ha.month)) 
hist(sqrt(LeafLitt2$Branches.Mg.ha.month))
hist(sqrt(LeafLitt2$Detritus.Mg.ha.month))
hist(sqrt(LeafLitt2$Flowers.Mg.ha.month))
hist(sqrt(LeafLitt2$Fruits.Mg.ha.month))
hist(sqrt(LeafLitt2$Total.Mg.ha.month))

Square-root seems to work well for Leaf, Detritus, and Total. Branches sorta.

Second test: Log transformation.

hist(log(LeafLitt2$Flowers.Mg.ha.month)) # Great
hist(log(LeafLitt2$Fruits.Mg.ha.month)) # Ok
hist(log(LeafLitt2$Branches.Mg.ha.month)) # Right skew.

Log works well for flowers and fruit. Branches needs something else.

hist(LeafLitt2$Branches.Mg.ha.month^(1/3)) 

Cube root best for Branches.

Note: Using the models without transformations does not yield radically different conclusions. However, we maintain the use of transformations to remain consistent with the literature and with reccommendations from biostatisticians.


Descriptive Results

Flowering Peaks

Flower litter, when does it fall with greater abundance? Let’s find out for each age group in each landscape. First for Hills. Let’s calculate means of Flower production by age Month using the aggregate function wrapped within a kbl function for style.

options(digits=7)
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = LomeRL, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()
Sample Flowers.Mg.ha.month
Aug 0.02000
Sept 0.00150
Oct 0.00250
Nov 0.01000
Dec 0.00025
Jan 0.00025
Feb 0.00100
Mar 0.00050
Apr 0.00075
May 0.00025
Jun 0.00200
Jul 0.00350
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = LomeDL2, FUN = mean, na.rm = TRUE), format = "html", 
table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Flowers.Mg.ha.month
Aug 0.0018333
Sept 0.0006667
Oct 0.0015000
Nov 0.0013333
Dec 0.0025000
Jan 0.0010000
Feb 0.0055000
Mar 0.0006667
Apr 0.0028333
May 0.0001667
Jun 0.0000000
Jul 0.0021667
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = LomeDL1, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Flowers.Mg.ha.month
Aug 0.0000000
Sept 0.0000000
Oct 0.0000000
Nov 0.0000000
Dec 0.0003333
Jan 0.0000000
Feb 0.0006667
Mar 0.0000000
Apr 0.0000000
May 0.0003333
Jun 0.0000000
Jul 0.0000000
ggplot(data=Lomerio, aes(x=Sample, y=Flowers.Mg.ha.month )) + # Flowers
  geom_bar(stat="identity", aes(fill = Succession), position="dodge") + 
  xlab("Month") +
  ylab(bquote('Flower Litter '(Mg~ha^-1~m^-1))) + 
  scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) + 
  guides(fill=guide_legend(title="Age Group")) +
  theme_bw(base_size = 12)

Next for the mountain landscape.

kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = MounRL, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Flowers.Mg.ha.month
Sept 0.0146667
Oct 0.0098333
Nov 0.0020000
Dec 0.0010000
Jan 0.0037500
Feb 0.0000000
Mar 0.0028333
Apr 0.0010000
May 0.0008333
Jun 0.0036667
Jul 0.0051667
Aug 0.0181667
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = MounDM2, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Flowers.Mg.ha.month
Sept 0.0026667
Oct 0.0035000
Nov 0.0036667
Dec 0.0035000
Jan 0.0010000
Feb 0.0018333
Mar 0.0038333
Apr 0.0005000
May 0.0075000
Jun 0.0063333
Jul 0.0023333
Aug 0.0015000
kbl(aggregate(Flowers.Mg.ha.month ~ Sample, data = MounDM1, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Flowers.Mg.ha.month
Sept 0.0010000
Oct 0.0001667
Nov 0.0003333
Dec 0.0018333
Jan 0.0091667
Feb 0.0008333
Mar 0.0016667
Apr 0.0000000
May 0.0005000
Jun 0.0005000
Jul 0.0021667
Aug 0.0008333
ggplot(data=Mountain, aes(x=Sample, y=Flowers.Mg.ha.month )) + # Flowers
  geom_bar(stat="identity", aes(fill = Succession), position="dodge") + 
  xlab("Month") +
  ylab(bquote('Flower Litter '(Mg~ha^-1~m^-1))) + 
  scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) + 
  guides(fill=guide_legend(title="Age Group")) +
  theme_bw(base_size = 10)

Fruiting Peaks

options(digits=5)
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = LomeRL, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()
Sample Fruits.Mg.ha.month
Aug 0.00400
Sept 0.02325
Oct 0.00925
Nov 0.02200
Dec 0.01150
Jan 0.02550
Feb 0.02700
Mar 0.03500
Apr 0.01650
May 0.01925
Jun 0.02450
Jul 0.02050
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = LomeDL2, FUN = mean, na.rm = TRUE), format = "html", 
table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Fruits.Mg.ha.month
Aug 0.00783
Sept 0.05217
Oct 0.02217
Nov 0.00967
Dec 0.01000
Jan 0.05083
Feb 0.03567
Mar 0.01717
Apr 0.02867
May 0.01883
Jun 0.00933
Jul 0.04167
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = LomeDL1, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Fruits.Mg.ha.month
Aug 0.00300
Sept 0.00167
Oct 0.00500
Nov 0.00333
Dec 0.01367
Jan 0.00667
Feb 0.00467
Mar 0.01433
Apr 0.02933
May 0.04300
Jun 0.00767
Jul 0.08167
ggplot(data=Lomerio, aes(x=Sample, y=Fruits.Mg.ha.month )) + # Fruits
  geom_bar(stat="identity", aes(fill = Succession), position="dodge") + 
  xlab("Month") +
  ylab(bquote('Fruit Litter '(Mg~ha^-1~m^-1))) + 
  scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) + 
  guides(fill=guide_legend(title="Age Group")) +
  theme_bw(base_size = 12)

Next for the mountain landscape.

kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = MounRL, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Fruits.Mg.ha.month
Sept 0.04417
Oct 0.04500
Nov 0.02167
Dec 0.01583
Jan 0.01225
Feb 0.01567
Mar 0.01833
Apr 0.02017
May 0.02583
Jun 0.01867
Jul 0.01383
Aug 0.04650
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = MounDM2, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Fruits.Mg.ha.month
Sept 0.06767
Oct 0.02100
Nov 0.01017
Dec 0.01200
Jan 0.01080
Feb 0.01483
Mar 0.01000
Apr 0.01317
May 0.02917
Jun 0.02350
Jul 0.01333
Aug 0.02900
kbl(aggregate(Fruits.Mg.ha.month ~ Sample, data = MounDM1, FUN = mean, na.rm = TRUE, na.action=na.exclude), format = "html", table.attr = "style='width:30%;'") %>% kable_classic()  
Sample Fruits.Mg.ha.month
Sept 0.00683
Oct 0.00283
Nov 0.00333
Dec 0.00133
Jan 0.00617
Feb 0.01167
Mar 0.01050
Apr 0.00567
May 0.00733
Jun 0.00933
Jul 0.00700
Aug 0.02450
ggplot(data=Mountain, aes(x=Sample, y=Fruits.Mg.ha.month )) + # Fruits
  geom_bar(stat="identity", aes(fill = Succession), position="dodge") + 
  xlab("Month") +
  ylab(bquote('Fruit Litter '(Mg~ha^-1~m^-1))) + 
  scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) + 
  guides(fill=guide_legend(title="Age Group")) +
  theme_bw(base_size = 10)

Reference Age Groups

Lomerio reference plots appear to have higher leaf litter and structural variable totals than Mountain reference plots. Let’s check. First we’ll create a new dataframe that holds only reference plots and we’ll call it ReferencePlots. Then, we will use the aggregate function to obtain means for each variable of interest, setting Age Groups (Succession) as the grouping variable.

options(digits=5)
ReferencePlots <- droplevels(subset(LeafLitt2, LeafLitt2$Succession == c("RM", "RL")))
kbl(aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month) ~ Succession,data = ReferencePlots, FUN = mean, na.action=na.exclude), format = "html", table.attr = "style='width:100%;'", row.names = FALSE) %>% 
  kable_classic()
Succession LeafLitter.Mg.ha.month Branches.Mg.ha.month Flowers.Mg.ha.month Fruits.Mg.ha.month Detritus.Mg.ha.month Total.Mg.ha.month
RL 0.75613 0.05987 0.00267 0.02425 0.18196 1.03412
RM 0.46803 0.05847 0.00494 0.02329 0.14229 0.70279
ReferenceStrucPlots <- droplevels(subset(Struc, Struc$Succession == c("RM", "RL")))
kbl(aggregate(cbind(TreesHA, BasalHA, Biomass, Richness, Simpson_1.D) ~ Landscape,
          data = ReferenceStrucPlots, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:100%;'", row.names = FALSE) %>% kable_classic()
Landscape TreesHA BasalHA Biomass Richness Simpson_1.D
Lomerio 614.00 25.35 252.05 87.500 0.9755
Mountain 438.67 13.00 117.67 43.667 0.9300

Sure enough, Lomerio (RL) has higher leaf, detritus, and total leaf litter than mountain (RM). RL also has higher tree density, basal area, biomass, richness, and a lil bit of diversity.

Rainfall Differences

How do plots differ in rainfall? We can take a look using the same combination of functions as before. First we’ll take a look at the Hill landscape plots.

kbl(aggregate(cbind(RainNow) ~ Plot, data = Lomerio, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:20%;'", row.names = FALSE) %>% kable_classic()
Plot RainNow
P10 308.03
P11 302.92
P12 308.61
P13 308.61
P14 324.03
P16 323.96
P2 313.99
P3 308.44
P4 296.90
P5 302.92
P6 313.42

Ok, in the Hillplots, rainfall dances around 300, a couple have 320 mm. Now for the mountain.

kbl(aggregate(cbind(RainNow) ~ Plot, data = Mountain, FUN = mean, na.rm = TRUE), format = "html", table.attr = "style='width:20%;'", row.names = FALSE) %>% kable_classic()
Plot RainNow
P1 256.25
P10 286.60
P11 288.45
P12 284.69
P13 280.93
P14 284.75
P15 265.91
P16 256.25
P17 265.91
P18 255.08
P2 256.25
P6 256.25
P7 256.25
P8 265.91
P9 265.91

In the mountain landscape, precipitation (mm) values dance around 260, a couple have 288.

So now, do the landscapes differ in precipitation? Let’s find out via anova.

anova(lm(data=LeafLitt2, RainNow ~ Paisaje))
## Analysis of Variance Table
## 
## Response: RainNow
##            Df  Sum Sq Mean Sq F value Pr(>F)   
## Paisaje     1  140235  140235    9.33 0.0024 **
## Residuals 370 5561807   15032                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Indeed, landscapes signficantly differ in the amount of rainfall they receive. Hill landscape likely receives more, let’s see.

kbl(summarySE(LeafLitt2, measurevar="RainNow", groupvars=c("Paisaje"), na.rm = TRUE)[,-4][,-5][,-7], format = "html", table.attr = "style='width:40%;'", row.names = FALSE) %>% kable_classic()
Paisaje N RainNow se
Lomerio 156 310.33 10.7280
Mountain 216 270.98 7.7356

Sunlight and Irradiation

Before anything, we need to get our Sunlight and Irradiation dataframes into formats that are easier to use. For this, we apply the melt function to create SolarMelt, a dataframe where we have a single column for Month, a single column for Municipio, and a single column for Sunlight.

SolarMelt <- melt(Solar,id.vars=c("Estacion", "Municipio"))
colnames(SolarMelt) <- c("Estacion", "Municipio", "Month", "Sunlight")

Now let’s visualize using ggplot. Do different Municipios in Caqueta receive wildly different hours of sunlight? Likely no.

ggplot(SolarMelt, aes(x=Month, y=Sunlight, group = Estacion, color = Estacion)) +
  geom_line(size=1, alpha=0.9, linetype=1) +
  ylab("Sunlight (hrs/day)") + 
  scale_colour_manual(values = c("#00AFBB", "#E7B800", "#FC4E07","#FFDB6D", "#C4961A", "#F4EDCA", 
                                 "#D16103", "#C3D7A4", "#52854C", "#4E84C4", "#293352")) + 
  theme_bw(base_size = 12) 

Great. It appears that all municipios share the same general sunlight pattern, so maybe we can simply use an average! Let’s extract the average sunlight per station and save it as SolarMeltAve.

SolarMeltAve <- droplevels(subset(SolarMelt, SolarMelt$Estacion == "Average"))

Next, we need to do some setup to create a dataframe that holds litter production and our climatic variables. First, we will use the select function from the dplyr package to specify the columns we wish to keep in our new dataframes, Lome3 and Moun3. Next, we want to calculate the average litter production and precipitation by month. These will then be Lome3Ave and Moun3Ave. After renaming columns with colnames, we merge with our sunlight dataframe using the merge function, creating the AveLom and AveMon dataframes.

Lome3 <- Lomerio %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation) 
Moun3 <- Mountain %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation)

Lome3Ave <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month, Precipitation) ~ Sample,
                      data = Lome3, FUN = mean, na.action=na.exclude)

Moun3Ave <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Detritus.Mg.ha.month, Total.Mg.ha.month, Precipitation) ~ Sample,
                      data = Moun3, FUN = mean, na.action=na.exclude)

colnames(Lome3Ave) <- c("Month", "LeafLitter.Mg.ha.month", "Branches.Mg.ha.month", "Detritus.Mg.ha.month", 
                        "Flowers.Mg.ha.month", "Fruits.Mg.ha.month", "Total.Mg.ha.month", "Precipitation")

colnames(Moun3Ave) <- c("Month", "LeafLitter.Mg.ha.month", "Branches.Mg.ha.month", "Detritus.Mg.ha.month", 
                        "Flowers.Mg.ha.month", "Fruits.Mg.ha.month", "Total.Mg.ha.month", "Precipitation")

AveLom <- merge(SolarMeltAve, Lome3Ave,  by="Month")
AveMon <- merge(SolarMeltAve, Moun3Ave, by="Month")

Next, we merge these dataframes with the Irad dataframe also through the merge function.

AveLomX <- merge(AveLom, Irad,  by="Month")
AveMonX <- merge(AveMon, Irad,  by="Month")

Now lets quickly visualize how precipitation and sunlight behave within the hill and mountain landscapes.

ggplot(AveLom, aes(x=Precipitation, y=Sunlight)) +
  geom_point(size=1, alpha=0.9, linetype=1) +
  geom_smooth(method=lm, colour = "black", fill="#23431F") +
  ylab("Sunlight (hrs/day)") + 
  theme_bw(base_size = 15) 

ggplot(AveMon, aes(x=Precipitation, y=Sunlight)) +
  geom_point(size=1, alpha=0.9, linetype=1) +
  geom_smooth(method=lm, colour = "black", fill="#67360D") +
  ylab("Sunlight (hrs/day)") + 
  theme_bw(base_size = 15) 

Good! Certainly look negatively correlated, but lets check using the cor.test function. We’ll run both, Pearson’s and Spearman’s Rank Correlation to check. We’ll also run the correaltions between precipitation and Irradiation.

cor.test(AveLomX$Precipitation, AveLomX$Sunlight, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  AveLomX$Precipitation and AveLomX$Sunlight
## t = -2.57, df = 10, p-value = 0.028
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.884504 -0.089213
## sample estimates:
##      cor 
## -0.63082
cor.test(AveLomX$Precipitation, AveLomX$Sunlight, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  AveLomX$Precipitation and AveLomX$Sunlight
## S = 460, p-value = 0.036
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.60702
cor.test(AveMonX$Precipitation, AveMonX$Sunlight, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  AveMonX$Precipitation and AveMonX$Sunlight
## t = -2.2, df = 10, p-value = 0.052
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8624520  0.0038045
## sample estimates:
##      cor 
## -0.57134
cor.test(AveMonX$Precipitation, AveMonX$Sunlight, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  AveMonX$Precipitation and AveMonX$Sunlight
## S = 453, p-value = 0.047
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.58246
cor.test(AveLomX$Precipitation, AveLomX$Irradiation, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  AveLomX$Precipitation and AveLomX$Irradiation
## t = -2.55, df = 10, p-value = 0.029
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.883188 -0.083241
## sample estimates:
##      cor 
## -0.62718
cor.test(AveLomX$Precipitation, AveLomX$Irradiation, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  AveLomX$Precipitation and AveLomX$Irradiation
## S = 452, p-value = 0.052
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.58042
cor.test(AveMonX$Precipitation, AveMonX$Irradiation, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  AveMonX$Precipitation and AveMonX$Irradiation
## t = -2.33, df = 10, p-value = 0.042
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.870921 -0.030231
## sample estimates:
##      cor 
## -0.59383
cor.test(AveMonX$Precipitation, AveMonX$Irradiation, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  AveMonX$Precipitation and AveMonX$Irradiation
## S = 454, p-value = 0.049
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.58741

Marvelous. Now, I need this to be separaed by habitat so I can illustrate it. First, we will use the select function from the dplyr package to specify the columns we wish to keep in our new dataframes, Lome4 and Moun4. We will then use the aggregate function to calculate means by Month (Sample) and Age Group (Succession), creating Lome4Ave and Moun4Ave. After renaming with colnames, we merge with SolarMeltAve so that we have values of sunlight in our new dataframes as well. We do this once more with the Irad dataframe to include irradiation. Note that this only works because all these dataframes are ordered in the same way.

Lome4 <- Lomerio %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation, Succession)
Moun4 <- Mountain %>% dplyr::select(Sample, LeafLitter.Mg.ha.month, Branches.Mg.ha.month, Detritus.Mg.ha.month, Flowers.Mg.ha.month, Fruits.Mg.ha.month, Total.Mg.ha.month, Precipitation, Succession)

Lome4Ave <- aggregate(Lome4[ , 2:8], by = list(Lome4$Sample, Lome4$Succession), FUN = mean)
Moun4Ave <- aggregate(Moun4[ , 2:8], by = list(Moun4$Sample, Moun4$Succession), FUN = mean, na.rm = TRUE)

colnames(Lome4Ave) <- c("Month", "Habitat", "LeafLitter.Mg.ha.month", 
                        "Branches.Mg.ha.month", "Detritus.Mg.ha.month", 
                        "Flowers.Mg.ha.month", "Fruits.Mg.ha.month", 
                        "Total.Mg.ha.month", "Precipitation")

colnames(Moun4Ave) <- c("Month",  "Habitat", "LeafLitter.Mg.ha.month", 
                        "Branches.Mg.ha.month", "Detritus.Mg.ha.month", 
                        "Flowers.Mg.ha.month", "Fruits.Mg.ha.month", 
                        "Total.Mg.ha.month", "Precipitation")

summaryAveLom <- merge(SolarMeltAve, Lome4Ave, by="Month")
summaryAveMon <- merge(SolarMeltAve, Moun4Ave, by="Month")

summaryAveLomX <- merge(summaryAveLom, Irad, by="Month") 
summaryAveMonX <- merge(summaryAveMon, Irad, by="Month") 

Next step is to set the Months into the proper format for figure making. Once relabeled, we can use the as.Date function to turn our variable into Date class. Then, we use the month function to make it so that our date class variable shows only the month. Finally, we rename the Age Groups into their final names.

levels(summaryAveLomX$Month) <- list("2017-08-01" = "Aug", "2017-09-01" =  "Sept", "2017-10-01" = "Oct", 
                                     "2017-11-01" = "Nov", "2017-12-01" = "Dec", "2018-01-01" = "Jan",
                                     "2018-02-01" = "Feb", "2018-03-01" = "Mar", "2018-04-01" = "Apr", 
                                     "2018-05-01" = "May", "2018-06-01" = "Jun","2018-07-01" = "Jul")

levels(summaryAveMonX$Month) <- list("2017-09-01" =  "Sept", "2017-10-01" = "Oct", 
                                     "2017-11-01" = "Nov", "2017-12-01" = "Dec", "2018-01-01" = "Jan",
                                     "2018-02-01" = "Feb", "2018-03-01" = "Mar", "2018-04-01" = "Apr", 
                                     "2018-05-01" = "May", "2018-06-01" = "Jun","2018-07-01" = "Jul",
                                     "2018-08-01" = "Aug")

summaryAveLomX$Month <- as.Date(summaryAveLomX$Month)
summaryAveLomX$Month <- month(summaryAveLomX$Month, label = TRUE)
summaryAveMonX$Month <- as.Date(summaryAveMonX$Month)
summaryAveMonX$Month <- month(summaryAveMonX$Month, label = TRUE)

levels(summaryAveLomX$Habitat) <- list(EH  = "DL1", IH = "DL2", OH = "RL")
levels(summaryAveMonX$Habitat) <- list(EM  = "DM1", IM = "DM2", OM = "RM")

We’re almost ready to graph! However, this complicated graph has three y-axes. I can only in R insert two, so we’ll have to make two versions of this graph that we at the end combine in PowerPoint.

Figure 4

This first version will include only the Precipitation as the secondary y-axis. We will create a graph using ggplot for each landscape and store them before combining them with a shared legend.

GLX1 <- ggplot(data=summaryAveLomX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
  geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) + 
  xlab("") +
  ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) + 
  geom_line(aes(Month, Irradiation*.003, group = Habitat), color = "GoldenRod") +
  geom_point(aes(Month, Irradiation*.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
  geom_boxplot(aes(Month, Precipitation*.002), color = "darkblue") + 
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
                     sec.axis = sec_axis(~./.002, name="Precipitation (mm)", breaks = seq(0,500,100))) + 
  scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) + 
  expand_limits(y = c(0, 1.1)) +   
  guides(fill=guide_legend(title="Age Group")) + 
  theme_bw(base_size = 12) + theme(legend.title = element_text(color = "black"))

GMX1 <- ggplot(data=summaryAveMonX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
  geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) + 
  xlab("Month") +
  ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) + 
  geom_line(aes(Month, Irradiation*.003, group = Habitat), color = "GoldenRod") +
  geom_point(aes(Month, Irradiation*.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
  geom_boxplot(aes(Month, Precipitation*.002), color = "darkblue") + 
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
                     sec.axis = sec_axis(~./0.002, name="Precipitation (mm)", breaks = seq(0,500,100))) +
  scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) + 
  expand_limits(y = c(0, 1.1)) +   
  guides(fill=guide_legend(title="")) +
  theme(legend.key.size = unit(0.2, 'cm')) +
  theme_bw(base_size = 12) 

combined1 <- GLX1 + GMX1 & theme(legend.position = "bottom") + 
  theme(legend.title = element_text(color = "black")) +  theme(legend.title.align=0.5) +
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))
combined1 + plot_layout(guides = "collect", nrow = 2) + plot_annotation(tag_levels = c('a'))

GLX2 <- ggplot(data=summaryAveLomX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
  geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) + 
  xlab("") +
  ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) + 
  geom_line(aes(Month, Irradiation*0.003, group = Habitat), color = "GoldenRod") +
  geom_point(aes(Month, Irradiation*0.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
  geom_boxplot(aes(Month, Precipitation*0.002), color = "darkblue") + 
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
                     sec.axis = sec_axis(~./0.003, name= bquote('Irradiation'~(W~m^-2)))) + 
  scale_fill_manual(values=c("#82C27B", "#43833C", "#23431F")) + 
  expand_limits(y = c(0, 1.1)) +   
  guides(fill=guide_legend(title="Age Group")) + 
  theme_bw(base_size = 12) + theme(legend.title = element_text(color = "black"))

GMX2 <- ggplot(data=summaryAveMonX, aes(x=Month, y=LeafLitter.Mg.ha.month )) +
  geom_bar(stat="identity", aes(fill = Habitat), position="dodge", alpha = 0.8) + 
  xlab("Month") +
  ylab(bquote('Leaf Litter '(Mg~ha^-1~mth^-1))) + 
  geom_line(aes(Month, Irradiation*0.003, group = Habitat), color = "GoldenRod") +
  geom_point(aes(Month, Irradiation*0.003, group = Habitat), color = "Black", shape = "☼", size = 3) +
  geom_boxplot(aes(Month, Precipitation*0.002), color = "darkblue") + 
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),
                     sec.axis = sec_axis(~./0.003, name= bquote('Irradiation'~(W~m^-2)))) + 
  scale_fill_manual(values=c("#FBBA84", "#CE803F", "#67360D")) + 
  expand_limits(y = c(0, 1.1)) +   
  guides(fill=guide_legend(title="")) +
  theme_bw(base_size = 12) + theme(legend.title = element_text(color = "black"))

combined2 <- GLX2 + GMX2 & theme(legend.position = "bottom") + 
  theme(legend.title = element_text(color = "black")) +  theme(legend.title.align=0.5) +
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))

combined2 + plot_layout(guides = "collect", nrow = 2) + plot_annotation(tag_levels = c('a'))

Statistical Analysis

Litter x Landscape

Does litterfall differ between lomerio and mountain? We can find out through Analysis of Variance (ANOVA) available through the anova function. We will use the appropriate transformations for each variable. For interpretation, typically a p-value below 0.05 is considered statistically significant, and a p-value between 0.05 and 0.1 as marginally significant.

anova(lm(data = LeafLitt2, sqrt(LeafLitter.Mg.ha.month) ~ Paisaje,na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(LeafLitter.Mg.ha.month)
##            Df Sum Sq Mean Sq F value Pr(>F)    
## Paisaje     1   1.29   1.292      47  3e-11 ***
## Residuals 367  10.09   0.027                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(data = LeafLitt2, Branches.Mg.ha.month^1/3 ~ Paisaje,na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: Branches.Mg.ha.month^1/3
##            Df Sum Sq  Mean Sq F value Pr(>F)   
## Paisaje     1 0.0012 0.001178    7.51 0.0064 **
## Residuals 367 0.0575 0.000157                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(data = LeafLitt2, log(Fruits.Mg.ha.month+0.05) ~ Paisaje,na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: log(Fruits.Mg.ha.month + 0.05)
##            Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje     1    0.2   0.156     1.7   0.19
## Residuals 367   33.7   0.092
anova(lm(data = LeafLitt2, log(Flowers.Mg.ha.month+0.05) ~ Paisaje,na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: log(Flowers.Mg.ha.month + 0.05)
##            Df Sum Sq Mean Sq F value Pr(>F)  
## Paisaje     1   0.05  0.0509    4.48  0.035 *
## Residuals 367   4.17  0.0114                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(data = LeafLitt2, sqrt(Detritus.Mg.ha.month) ~ Paisaje,na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(Detritus.Mg.ha.month)
##            Df Sum Sq Mean Sq F value Pr(>F)
## Paisaje     1   0.00 0.00017    0.01   0.91
## Residuals 367   4.82 0.01314
anova(lm(data = LeafLitt2, sqrt(Total.Mg.ha.month) ~ Paisaje,na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(Total.Mg.ha.month)
##            Df Sum Sq Mean Sq F value  Pr(>F)    
## Paisaje     1   0.84   0.842    23.3 2.1e-06 ***
## Residuals 367  13.28   0.036                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Figure 2

Lets make a stacked barplot to visualize how leaf litter is composed in each landscape. To do this, first we need to make a new dataframe, LomerioStack1, that will hold only the columsn that we are interested in. For this, we can use the -which function to select the columns we DONT want by name. Yeah, we could have just used the which function and select the ones we want, but, we’re still learning!

LomerioStack1 <- Lomerio[ , -which(names(Lomerio) %in% c("Record","Paisaje", "Locality", 
                 "Successional.State", "Plot", 
                 "LeafLitter.g", "LeafLitter.Mg.ha.month2",
                 "Branches.g", "Detritus.g",
                 "Flowers.g", "Fruits.g", "Total.g",
                 "Total.Mg.ha.month", "Precipitation","Temp.Ave","Temp.Max", 
                 "PlotSample", "CodeX",  "CodeZ","CodeY", "RainNow",  
                 "Rain1",  "Rain2",  "Rain3", "TempNow", 
                 "Temp1", "Temp2", "Temp3"))]

Next we use the melt function from the reshape package to change the format of our table so that we have a new column that holds all litter types. We’ll call this new dataframe LomerioStack2. We will then use colnames to rename the columns into what we want. Next, we use the summarise function from the dplyr package to calculate mean values of litter production for each specific grouping. Why mean and not sum? Well, because our landscapes do not have an equal number of plots, and we’re more interested in the proportion of littertype. We’ll throw those mean values into a new dataframe, LomerioStack3.

LomerioStack2 <- melt(LomerioStack1, id.vars=c("Sample", "Succession"))
colnames(LomerioStack2) <- c("Month", "Habitat", "LitterType", "Value")

LomerioStack3 <- as.data.frame(LomerioStack2 %>% 
                                 group_by(LitterType, Habitat) %>% 
                                 dplyr::summarise(Value = mean(Value)))

We’ll do the exact same for the mountain plots, with the exception that halfway through we’ll use the na.omit function to remove a few NA values.

MountainStack1 <- Mountain[ , -which(names(Mountain) %in% c("Record","Paisaje", "Locality", "Successional.State", "Plot", "LeafLitter.g", "LeafLitter.Mg.ha.month2",
"Branches.g", "Detritus.g", "Flowers.g", "Fruits.g", "Total.g",
"Total.Mg.ha.month", "Precipitation","Temp.Ave","Temp.Max", "PlotSample",
"CodeX",  "CodeZ","CodeY", "RainNow",  "Rain1",  "Rain2",  "Rain3", "TempNow", 
"Temp1", "Temp2", "Temp3"))]

MountainStack2 <- melt(MountainStack1, id.vars=c("Sample", "Succession"))
colnames(MountainStack2) <- c("Month", "Habitat", "LitterType", "Value")

MountainStack2.5 <- na.omit(MountainStack2)

MountainStack3 <- as.data.frame(MountainStack2.5 %>% 
                                  group_by(LitterType, Habitat) %>% 
                                  dplyr::summarise(Value = mean(Value)))

Time to graph! We’ll use ggplot to create this stacked barplots, but first, we will rename the landscapes to their english names and bind them together by row using the rbind function, creating the dataframe Megastack. We’ll also rename the Age Groups into the final names the publication needs. Finally, with colnames, we rename the columns.

LomerioStack3$Landscape <- "Hill"
MountainStack3$Landscape <- "Mountain"

Megastack <- rbind(LomerioStack3, MountainStack3)

levels(Megastack$Habitat) <- list(EH  = "DL1", IH = "DL2", OH = "RL",
                                  EM  = "DM1", IM = "DM2", OM = "RM")
colnames(Megastack) <- c("LitterType", "Age Group", "Value", "Landscape")
ggplot(Megastack, aes(fill=LitterType, y=Value, x=`Age Group`)) + 
  geom_bar(position="stack", stat="identity", color = "black", alpha = 0.7) + 
  ylab(bquote('Mean Total Litter '(Mg~ha^-1~mth^-1))) + 
  ylim(0, 1) + 
  labs(fill = "") + 
  scale_fill_brewer(palette = "Accent", labels=c("Leaf", "Branch", "Detritus", 
                                                 "Flower", "Fruit"),
                    guide = guide_legend(direction = "horizontal",
                      title.position = "left")) + 
  theme_bw(base_size = 13) + theme(legend.title = element_text(color = "black")) + 
  theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid')) +
  theme(legend.position = "bottom")  + 
  theme(legend.title.align=0.5) +
  theme(legend.key.size = unit(0.4, 'cm')) +
  facet_wrap(~Landscape, scales = "free") 

Litter x Precipitation

Do the habitats differ in their relationship between leaf litter and precipitation? We can find out by running an ANCOVA with Precipitation and Habitat as interactive terms in an anova function.

anova(lm(data = summaryAveLom, sqrt(LeafLitter.Mg.ha.month) ~ Precipitation * Habitat))
## Analysis of Variance Table
## 
## Response: sqrt(LeafLitter.Mg.ha.month)
##                       Df Sum Sq Mean Sq F value  Pr(>F)    
## Precipitation          1 0.0344  0.0344    3.94   0.056 .  
## Habitat                2 0.2788  0.1394   15.98 1.9e-05 ***
## Precipitation:Habitat  2 0.0030  0.0015    0.17   0.841    
## Residuals             30 0.2616  0.0087                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(data = summaryAveMon, sqrt(LeafLitter.Mg.ha.month) ~ Precipitation * Habitat))
## Analysis of Variance Table
## 
## Response: sqrt(LeafLitter.Mg.ha.month)
##                       Df Sum Sq Mean Sq F value  Pr(>F)    
## Precipitation          1 0.1527  0.1527   24.23 2.9e-05 ***
## Habitat                2 0.0789  0.0395    6.26  0.0053 ** 
## Precipitation:Habitat  2 0.0025  0.0012    0.20  0.8238    
## Residuals             30 0.1890  0.0063                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

They do not, so we will be pooling the habitats within landscape together.

Here we’re going to include interactions between Precipitation and Month as well as Habitat and Month.

Hill

anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: sqrt(LeafLitter.Mg.ha.month)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2  0.923   0.461   25.43 6.8e-10 ***
## RainNow              1  0.126   0.126    6.94  0.0096 ** 
## Sample              11  0.936   0.085    4.69 6.4e-06 ***
## Succession:RainNow   2  0.008   0.004    0.21  0.8076    
## Succession:Sample   22  0.351   0.016    0.88  0.6221    
## Residuals          117  2.122   0.018                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession + RainNow + Sample + Succession + Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: sqrt(LeafLitter.Mg.ha.month)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession   2  0.923   0.461   26.22 2.1e-10 ***
## RainNow      1  0.126   0.126    7.15  0.0084 ** 
## Sample      11  0.936   0.085    4.84 2.6e-06 ***
## Residuals  141  2.481   0.018                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: Branches.Mg.ha.month^(1/3)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2  0.220  0.1102   16.70 4.2e-07 ***
## RainNow              1  0.013  0.0131    1.99 0.16126    
## Sample              11  0.270  0.0245    3.72 0.00014 ***
## Succession:RainNow   2  0.002  0.0012    0.17 0.83982    
## Succession:Sample   22  0.124  0.0056    0.85 0.65681    
## Residuals          117  0.772  0.0066                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession + RainNow + Sample + Succession + Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: Branches.Mg.ha.month^(1/3)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession   2  0.220  0.1102   17.30 1.9e-07 ***
## RainNow      1  0.013  0.0131    2.06    0.15    
## Sample      11  0.270  0.0245    3.86 7.2e-05 ***
## Residuals  141  0.898  0.0064                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Flowers.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: log(Flowers.Mg.ha.month + 0.05)
##                     Df Sum Sq Mean Sq F value Pr(>F)   
## Succession           2  0.068  0.0338    5.61 0.0047 **
## RainNow              1  0.001  0.0005    0.08 0.7724   
## Sample              11  0.134  0.0122    2.02 0.0320 * 
## Succession:RainNow   2  0.028  0.0138    2.29 0.1060   
## Succession:Sample   22  0.247  0.0112    1.86 0.0185 * 
## Residuals          117  0.706  0.0060                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Flowers.Mg.ha.month+0.05) ~ Succession + RainNow + Sample + Succession + Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: log(Flowers.Mg.ha.month + 0.05)
##             Df Sum Sq Mean Sq F value Pr(>F)   
## Succession   2  0.068  0.0338    4.87 0.0091 **
## RainNow      1  0.001  0.0005    0.07 0.7876   
## Sample      11  0.134  0.0122    1.76 0.0675 . 
## Residuals  141  0.980  0.0070                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: log(Fruits.Mg.ha.month + 0.05)
##                     Df Sum Sq Mean Sq F value Pr(>F)  
## Succession           2   0.33  0.1662    1.61   0.21  
## RainNow              1   0.09  0.0869    0.84   0.36  
## Sample              11   1.94  0.1765    1.71   0.08 .
## Succession:RainNow   2   0.23  0.1171    1.13   0.33  
## Succession:Sample   22   1.64  0.0746    0.72   0.81  
## Residuals          117  12.11  0.1035                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession+ RainNow + Sample + Succession + Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: log(Fruits.Mg.ha.month + 0.05)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Succession   2   0.33  0.1662    1.68  0.191  
## RainNow      1   0.09  0.0869    0.88  0.351  
## Sample      11   1.94  0.1765    1.78  0.063 .
## Residuals  141  13.98  0.0992                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Detritus.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: sqrt(Detritus.Mg.ha.month)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2  0.311  0.1555   17.57 2.1e-07 ***
## RainNow              1  0.104  0.1037   11.71 0.00086 ***
## Sample              11  0.131  0.0119    1.34 0.20902    
## Succession:RainNow   2  0.052  0.0260    2.94 0.05676 .  
## Succession:Sample   22  0.196  0.0089    1.01 0.46200    
## Residuals          117  1.035  0.0089                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Detritus has marginally sig interaction Succession:RainNow.

anova(lm(sqrt(Total.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: sqrt(Total.Mg.ha.month)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2  1.382   0.691   31.71 9.9e-12 ***
## RainNow              1  0.236   0.236   10.84 0.00131 ** 
## Sample              11  0.892   0.081    3.72 0.00014 ***
## Succession:RainNow   2  0.027   0.013    0.61 0.54305    
## Succession:Sample   22  0.451   0.020    0.94 0.54392    
## Residuals          117  2.549   0.022                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Total.Mg.ha.month) ~ Succession+ RainNow + Sample + Succession + Sample, data = Lomerio))
## Analysis of Variance Table
## 
## Response: sqrt(Total.Mg.ha.month)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession   2  1.382   0.691   32.19 3.1e-12 ***
## RainNow      1  0.236   0.236   11.00  0.0012 ** 
## Sample      11  0.892   0.081    3.78 9.3e-05 ***
## Residuals  141  3.027   0.021                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Because only Detritus had a marginally significant interaction, for detritus only we will separate by Age Group and evaluate the relationship between rain and detritus for each group. First, we use droplevels and subset to create these groups. Then, we run the anovas again.

LomeDL1 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL1"))
LomeDL2 <- droplevels(subset(Lomerio, Lomerio$Succession == "DL2"))
LomeRL <- droplevels(subset(Lomerio, Lomerio$Succession == "RL"))
anova(lm(sqrt(Detritus.Mg.ha.month) ~ RainNow + Sample, data = LomeDL1)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus.Mg.ha.month)
##           Df Sum Sq Mean Sq F value Pr(>F)
## RainNow    1 0.0033 0.00334    0.51   0.48
## Sample    11 0.0847 0.00770    1.17   0.36
## Residuals 23 0.1517 0.00659
anova(lm(sqrt(Detritus.Mg.ha.month) ~ RainNow + Sample, data = LomeDL2)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus.Mg.ha.month)
##           Df Sum Sq Mean Sq F value Pr(>F)   
## RainNow    1  0.097  0.0973   10.33 0.0021 **
## Sample    11  0.129  0.0117    1.25 0.2775   
## Residuals 59  0.555  0.0094                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Detritus.Mg.ha.month) ~ RainNow + Sample, data = LomeRL)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus.Mg.ha.month)
##           Df Sum Sq Mean Sq F value Pr(>F)  
## RainNow    1  0.057  0.0569    6.06  0.019 *
## Sample    11  0.111  0.0101    1.08  0.407  
## Residuals 35  0.329  0.0094                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Next we use the cor.test function to obtain Pearson’s and Spearman’s Rank coefficients.

cor.test(sqrt(Lomerio$LeafLitter.Mg.ha.month), Lomerio$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(Lomerio$LeafLitter.Mg.ha.month) and Lomerio$RainNow
## t = -2.07, df = 154, p-value = 0.041
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3132760 -0.0072663
## sample estimates:
##      cor 
## -0.16422
cor.test(Lomerio$LeafLitter.Mg.ha.month, Lomerio$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Lomerio$LeafLitter.Mg.ha.month and Lomerio$RainNow
## S = 727198, p-value = 0.063
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.14934
cor.test(Lomerio$Branches.Mg.ha.month^(1/3), Lomerio$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  Lomerio$Branches.Mg.ha.month^(1/3) and Lomerio$RainNow
## t = -1.16, df = 154, p-value = 0.25
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.246880  0.064729
## sample estimates:
##       cor 
## -0.093361
cor.test(Lomerio$Branches.Mg.ha.month, Lomerio$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Lomerio$Branches.Mg.ha.month and Lomerio$RainNow
## S = 703803, p-value = 0.16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.11236
cor.test(log(Lomerio$Flowers.Mg.ha.month+0.05), Lomerio$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  log(Lomerio$Flowers.Mg.ha.month + 0.05) and Lomerio$RainNow
## t = -0.237, df = 154, p-value = 0.81
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.17569  0.13848
## sample estimates:
##       cor 
## -0.019076
cor.test(Lomerio$Flowers.Mg.ha.month, Lomerio$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Lomerio$Flowers.Mg.ha.month and Lomerio$RainNow
## S = 660023, p-value = 0.59
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.043168
cor.test(log(Lomerio$Fruits.Mg.ha.month+0.05), Lomerio$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  log(Lomerio$Fruits.Mg.ha.month + 0.05) and Lomerio$RainNow
## t = -0.892, df = 154, p-value = 0.37
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.226265  0.086442
## sample estimates:
##       cor 
## -0.071672
cor.test(Lomerio$Fruits.Mg.ha.month, Lomerio$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Lomerio$Fruits.Mg.ha.month and Lomerio$RainNow
## S = 682877, p-value = 0.33
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.079289
cor.test(sqrt(Lomerio$Total.Mg.ha.month), Lomerio$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(Lomerio$Total.Mg.ha.month) and Lomerio$RainNow
## t = -2.57, df = 154, p-value = 0.011
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.348510 -0.046806
## sample estimates:
##      cor 
## -0.20246
cor.test(Lomerio$Total.Mg.ha.month, Lomerio$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Lomerio$Total.Mg.ha.month and Lomerio$RainNow
## S = 742872, p-value = 0.03
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.17411
cor.test(sqrt(LomeDL1$Detritus.Mg.ha.month), LomeDL1$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(LomeDL1$Detritus.Mg.ha.month) and LomeDL1$RainNow
## t = 0.693, df = 34, p-value = 0.49
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.21905  0.42985
## sample estimates:
##     cor 
## 0.11797
cor.test(LomeDL1$Detritus.Mg.ha.month, LomeDL1$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  LomeDL1$Detritus.Mg.ha.month and LomeDL1$RainNow
## S = 7450, p-value = 0.81
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.041189
cor.test(sqrt(LomeDL2$Detritus.Mg.ha.month), LomeDL2$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(LomeDL2$Detritus.Mg.ha.month) and LomeDL2$RainNow
## t = -3.15, df = 70, p-value = 0.0024
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.54026 -0.13185
## sample estimates:
##      cor 
## -0.35274
cor.test(LomeDL2$Detritus.Mg.ha.month, LomeDL2$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  LomeDL2$Detritus.Mg.ha.month and LomeDL2$RainNow
## S = 82770, p-value = 0.0045
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## -0.3308
cor.test(sqrt(LomeRL$Detritus.Mg.ha.month), LomeRL$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(LomeRL$Detritus.Mg.ha.month) and LomeRL$RainNow
## t = -2.44, df = 46, p-value = 0.019
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.56804 -0.06021
## sample estimates:
##      cor 
## -0.33855
cor.test(LomeRL$Detritus.Mg.ha.month, LomeRL$RainNow, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  LomeRL$Detritus.Mg.ha.month and LomeRL$RainNow
## S = 25976, p-value = 0.0038
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.40988

Mountain

We’ll repeat the process for the mountain landscape.

anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(LeafLitter.Mg.ha.month)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2  0.531   0.266   15.03 9.5e-07 ***
## RainNow              1  0.781   0.781   44.17 3.7e-10 ***
## Sample              11  0.881   0.080    4.53 5.0e-06 ***
## Succession:RainNow   2  0.006   0.003    0.17    0.84    
## Succession:Sample   22  0.349   0.016    0.90    0.60    
## Residuals          174  3.075   0.018                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(LeafLitter.Mg.ha.month) ~ Succession + RainNow + Sample + Succession + Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(LeafLitter.Mg.ha.month)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession   2   0.53   0.266   15.34 6.4e-07 ***
## RainNow      1   0.78   0.781   45.07 2.0e-10 ***
## Sample      11   0.88   0.080    4.62 2.9e-06 ***
## Residuals  198   3.43   0.017                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession* RainNow + Sample + Succession * Sample, data =  Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: Branches.Mg.ha.month^(1/3)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2  0.328  0.1638   21.34 5.1e-09 ***
## RainNow              1  0.001  0.0015    0.19    0.66    
## Sample              11  0.313  0.0285    3.71 9.2e-05 ***
## Succession:RainNow   2  0.006  0.0031    0.40    0.67    
## Succession:Sample   22  0.146  0.0066    0.86    0.64    
## Residuals          174  1.335  0.0077                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Branches.Mg.ha.month^(1/3) ~ Succession+ RainNow + Sample + Succession + Sample, data =  Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: Branches.Mg.ha.month^(1/3)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession   2  0.328  0.1638   21.80 2.8e-09 ***
## RainNow      1  0.001  0.0015    0.20    0.66    
## Sample      11  0.313  0.0285    3.79 6.0e-05 ***
## Residuals  198  1.487  0.0075                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Flowers.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: log(Flowers.Mg.ha.month + 0.05)
##                     Df Sum Sq Mean Sq F value Pr(>F)   
## Succession           2  0.119  0.0593    4.94 0.0082 **
## RainNow              1  0.001  0.0009    0.07 0.7849   
## Sample              11  0.218  0.0198    1.65 0.0897 . 
## Succession:RainNow   2  0.060  0.0298    2.48 0.0871 . 
## Succession:Sample   22  0.501  0.0228    1.89 0.0125 * 
## Residuals          174  2.091  0.0120                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: log(Fruits.Mg.ha.month + 0.05)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2   1.80   0.902   13.22 4.5e-06 ***
## RainNow              1   0.20   0.202    2.96  0.0871 .  
## Sample              11   2.26   0.205    3.01  0.0011 ** 
## Succession:RainNow   2   0.21   0.104    1.53  0.2195    
## Succession:Sample   22   1.06   0.048    0.70  0.8313    
## Residuals          174  11.87   0.068                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Fruits.Mg.ha.month+0.05) ~ Succession+ RainNow + Sample + Succession + Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: log(Fruits.Mg.ha.month + 0.05)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession   2   1.80   0.902   13.60 2.9e-06 ***
## RainNow      1   0.20   0.202    3.04 0.08257 .  
## Sample      11   2.26   0.205    3.09 0.00074 ***
## Residuals  198  13.14   0.066                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Detritus.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(Detritus.Mg.ha.month)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2  0.533  0.2663   25.48   2e-10 ***
## RainNow              1  0.013  0.0129    1.24 0.26728    
## Sample              11  0.362  0.0329    3.15 0.00066 ***
## Succession:RainNow   2  0.010  0.0052    0.50 0.60784    
## Succession:Sample   22  0.255  0.0116    1.11 0.34213    
## Residuals          174  1.819  0.0105                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Detritus.Mg.ha.month) ~ Succession+ RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(Detritus.Mg.ha.month)
##                    Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession          2  0.533  0.2663   25.15 2.5e-10 ***
## RainNow             1  0.013  0.0129    1.22 0.27033    
## Sample             11  0.362  0.0329    3.11 0.00076 ***
## Succession:Sample  22  0.221  0.0100    0.95 0.53420    
## Residuals         176  1.863  0.0106                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Total.Mg.ha.month) ~ Succession* RainNow + Sample + Succession * Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(Total.Mg.ha.month)
##                     Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession           2   1.33   0.665   29.43 9.8e-12 ***
## RainNow              1   0.59   0.586   25.94 9.1e-07 ***
## Sample              11   1.33   0.121    5.34 2.8e-07 ***
## Succession:RainNow   2   0.00   0.001    0.06    0.94    
## Succession:Sample   22   0.56   0.025    1.12    0.33    
## Residuals          174   3.93   0.023                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Total.Mg.ha.month) ~ Succession+ RainNow + Sample + Succession + Sample, data = Mountain, na.action=na.exclude))
## Analysis of Variance Table
## 
## Response: sqrt(Total.Mg.ha.month)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession   2   1.33   0.665   29.32 7.1e-12 ***
## RainNow      1   0.59   0.586   25.84 8.6e-07 ***
## Sample      11   1.33   0.121    5.32 2.3e-07 ***
## Residuals  198   4.49   0.023                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.test(sqrt(Mountain$LeafLitter.Mg.ha.month), Mountain$RainNow,  method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(Mountain$LeafLitter.Mg.ha.month) and Mountain$RainNow
## t = -5.76, df = 211, p-value = 3e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.47906 -0.24611
## sample estimates:
##      cor 
## -0.36835
cor.test(Mountain$LeafLitter.Mg.ha.month, Mountain$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Mountain$LeafLitter.Mg.ha.month and Mountain$RainNow
## S = 2164663, p-value = 2.6e-07
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.34404
cor.test(Mountain$Branches.Mg.ha.month^(1/3), Mountain$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  Mountain$Branches.Mg.ha.month^(1/3) and Mountain$RainNow
## t = -0.335, df = 211, p-value = 0.74
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.15700  0.11172
## sample estimates:
##       cor 
## -0.023054
cor.test(Mountain$Branches.Mg.ha.month, Mountain$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Mountain$Branches.Mg.ha.month and Mountain$RainNow
## S = 1600678, p-value = 0.93
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.0061382
cor.test(log(Mountain$Flowers.Mg.ha.month+0.05), Mountain$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  log(Mountain$Flowers.Mg.ha.month + 0.05) and Mountain$RainNow
## t = 0.315, df = 211, p-value = 0.75
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.11311  0.15563
## sample estimates:
##      cor 
## 0.021655
cor.test(Mountain$Flowers.Mg.ha.month, Mountain$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Mountain$Flowers.Mg.ha.month and Mountain$RainNow
## S = 1674001, p-value = 0.57
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.039388
cor.test(log(Mountain$Fruits.Mg.ha.month+0.05), Mountain$RainNow, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  log(Mountain$Fruits.Mg.ha.month + 0.05) and Mountain$RainNow
## t = 1.65, df = 211, p-value = 0.1
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.021708  0.243780
## sample estimates:
##     cor 
## 0.11305
cor.test(Mountain$Fruits.Mg.ha.month, Mountain$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Mountain$Fruits.Mg.ha.month and Mountain$RainNow
## S = 1456977, p-value = 0.17
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.095362
cor.test(sqrt(Mountain$Detritus.Mg.ha.month), Mountain$RainNow,  method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(Mountain$Detritus.Mg.ha.month) and Mountain$RainNow
## t = -0.908, df = 211, p-value = 0.37
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.195169  0.072667
## sample estimates:
##       cor 
## -0.062374
cor.test(Mountain$Detritus.Mg.ha.month, Mountain$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Mountain$Detritus.Mg.ha.month and Mountain$RainNow
## S = 1700554, p-value = 0.42
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.055875
cor.test(sqrt(Mountain$Total.Mg.ha.month), Mountain$RainNow,  method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  sqrt(Mountain$Total.Mg.ha.month) and Mountain$RainNow
## t = -4.08, df = 211, p-value = 6.3e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.39077 -0.14126
## sample estimates:
##      cor 
## -0.27055
cor.test(Mountain$Total.Mg.ha.month, Mountain$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Mountain$Total.Mg.ha.month and Mountain$RainNow
## S = 2018308, p-value = 0.00019
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.25317

We had a marginally significant interaction with flowers, so lets separate.

MounDM1 <- droplevels(subset(Mountain, Mountain$Succession == "DM1"))
MounDM2 <- droplevels(subset(Mountain, Mountain$Succession == "DM2"))
MounRM <- droplevels(subset(Mountain, Mountain$Succession == "RM"))

anova(lm(log(Flowers.Mg.ha.month+0.05) ~ RainNow + Sample, data = MounDM1)) # 0.08966
## Analysis of Variance Table
## 
## Response: log(Flowers.Mg.ha.month + 0.05)
##           Df Sum Sq Mean Sq F value Pr(>F)  
## RainNow    1  0.016 0.01597    2.98   0.09 .
## Sample    11  0.089 0.00808    1.51   0.15  
## Residuals 59  0.317 0.00536                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Flowers.Mg.ha.month+0.05) ~ RainNow + Sample, data = MounDM2)) # 0.06538
## Analysis of Variance Table
## 
## Response: log(Flowers.Mg.ha.month + 0.05)
##           Df Sum Sq Mean Sq F value Pr(>F)  
## RainNow    1  0.034  0.0342    3.53  0.065 .
## Sample    11  0.093  0.0085    0.88  0.568  
## Residuals 58  0.562  0.0097                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(log(Flowers.Mg.ha.month+0.05) ~ RainNow + Sample, data = MounRM)) # 0.98921
## Analysis of Variance Table
## 
## Response: log(Flowers.Mg.ha.month + 0.05)
##           Df Sum Sq Mean Sq F value Pr(>F)  
## RainNow    1  0.000  0.0000    0.00  0.989  
## Sample    11  0.546  0.0497    2.34  0.019 *
## Residuals 57  1.212  0.0213                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.test(MounDM1$Flowers.Mg.ha.month, MounDM1$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  MounDM1$Flowers.Mg.ha.month and MounDM1$RainNow
## S = 72221, p-value = 0.18
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.16119
cor.test(MounDM2$Flowers.Mg.ha.month, MounDM2$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  MounDM2$Flowers.Mg.ha.month and MounDM2$RainNow
## S = 59900, p-value = 0.97
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0043576
cor.test(MounRM$Flowers.Mg.ha.month, MounRM$RainNow,  method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  MounRM$Flowers.Mg.ha.month and MounRM$RainNow
## S = 56523, p-value = 0.93
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.011066

Litter x Structure

Do the litter traits correlate with the tree community variables? Before we’re able to answer this question, we need to calculate total litter and precipitation by our CodeZ variable. We can accomplish this with the aggregate function. We’ll then merge our dataframes together and rename using colnames.

DFX1 <- aggregate(cbind(LeafLitter.Mg.ha.month, Branches.Mg.ha.month, 
                        Flowers.Mg.ha.month, Fruits.Mg.ha.month, 
                        Detritus.Mg.ha.month, Total.Mg.ha.month) ~ CodeZ, 
                  data = LeafLitt2, FUN = sum, na.action=na.exclude)

DFX2 <- aggregate(cbind(RainNow, Rain1, Rain2, Rain3) ~ CodeZ,
                  data = LeafLitt2, FUN = mean, na.action=na.exclude)

DFX3 <- merge(DFX1, DFX2, by="CodeZ")

colnames(DFX3) <- c("CodeZ", "Leaf","Branch","Flower","Fruit","Detritus","Total",
                    "RainNow","Rain1","Rain2","Rain3")

Next we merge our new dataframe, DFX3, with the Struc dataframe that contains basal area, stem density, species density, biomass, and species diversity. We then subset to create a dataframe for Hill and one for Mountain. We’ll also subset again to create dataframes for specific age groups within landscape.

Struc <- merge(Struc, DFX3, by="CodeZ") # After merging we split by landscape.
StrucLome <- droplevels(subset(Struc, Struc$Landscape == "Lomerio"))
StrucMon <- droplevels(subset(Struc, Struc$Landscape == "Mountain"))

StrucDL1 <- droplevels(subset(Struc, Struc$Succession == "DL1")) # Lomerio
StrucDL2 <- droplevels(subset(Struc, Struc$Succession == "DL2"))
StrucRL <- droplevels(subset(Struc, Struc$Succession == "RL"))
StrucDM1 <- droplevels(subset(Struc, Struc$Succession == "DM1")) # Mountain
StrucDM2 <- droplevels(subset(Struc, Struc$Succession == "DM2"))
StrucRM <- droplevels(subset(Struc, Struc$Succession == "RM"))

Now we can look at correlations. We will use the non-parametric Spearman’s Rank Correlation for this.

Leaf Litter

First, we’ll check out the correlations without separating by Age Group. We’ll start with leaf litter.

cor.test(Struc$Leaf, Struc$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Leaf and Struc$TreesHA
## S = 4735, p-value = 0.81
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.045372
cor.test(Struc$Leaf, Struc$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Leaf and Struc$BasalHA
## S = 3743, p-value = 0.18
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.24531
cor.test(Struc$Leaf, Struc$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Leaf and Struc$Biomass
## S = 3718, p-value = 0.17
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.25043
cor.test(Struc$Leaf, Struc$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Leaf and Struc$Richness
## S = 3615, p-value = 0.14
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.27127
cor.test(Struc$Leaf, Struc$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Leaf and Struc$Simpson_1.D
## S = 3362, p-value = 0.077
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.32226

Next we will use ANCOVAs to determine if landscapes differ in their relationship between leaf litter and the structural variables.

anova(lm(sqrt(Leaf) ~ TreesHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Leaf)
##                   Df Sum Sq Mean Sq F value  Pr(>F)    
## TreesHA            1  0.083   0.083    1.68  0.2062    
## Landscape          1  1.730   1.730   35.04 2.6e-06 ***
## TreesHA:Landscape  1  0.590   0.590   11.95  0.0018 ** 
## Residuals         27  1.333   0.049                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Leaf) ~ BasalHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Leaf)
##                   Df Sum Sq Mean Sq F value Pr(>F)    
## BasalHA            1  0.226   0.226    4.01 0.0554 .  
## Landscape          1  1.501   1.501   26.59  2e-05 ***
## BasalHA:Landscape  1  0.486   0.486    8.61 0.0068 ** 
## Residuals         27  1.524   0.056                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Leaf) ~ Biomass * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Leaf)
##                   Df Sum Sq Mean Sq F value  Pr(>F)    
## Biomass            1  0.269   0.269    4.59   0.041 *  
## Landscape          1  1.460   1.460   24.97 3.1e-05 ***
## Biomass:Landscape  1  0.429   0.429    7.35   0.012 *  
## Residuals         27  1.578   0.058                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Leaf) ~ Richness * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Leaf)
##                    Df Sum Sq Mean Sq F value Pr(>F)    
## Richness            1  0.315   0.315    5.16  0.031 *  
## Landscape           1  1.412   1.412   23.16  5e-05 ***
## Richness:Landscape  1  0.362   0.362    5.94  0.022 *  
## Residuals          27  1.647   0.061                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Leaf) ~ Simpson_1.D * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Leaf)
##                       Df Sum Sq Mean Sq F value  Pr(>F)    
## Simpson_1.D            1  0.231   0.231    2.77 0.10735    
## Landscape              1  1.234   1.234   14.82 0.00066 ***
## Simpson_1.D:Landscape  1  0.021   0.021    0.25 0.62019    
## Residuals             27  2.250   0.083                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There was a significant interaction in all but diversity. Next we separate by landscape.

cor.test(StrucLome$Leaf, StrucLome$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Leaf and StrucLome$TreesHA
## S = 62, p-value = 0.00078
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.82967
cor.test(StrucLome$Leaf, StrucLome$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Leaf and StrucLome$Biomass
## S = 68, p-value = 0.0012
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.81319
cor.test(StrucLome$Leaf, StrucLome$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Leaf and StrucLome$Richness
## S = 59.6, p-value = 0.00037
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.83631
cor.test(StrucLome$Leaf, StrucLome$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Leaf and StrucLome$Simpson_1.D
## S = 78, p-value = 0.0023
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.78571
cor.test(StrucLome$Leaf, StrucLome$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Leaf and StrucLome$BasalHA
## S = 62, p-value = 0.00078
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.82967
cor.test(StrucMon$Leaf, StrucMon$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Leaf and StrucMon$TreesHA
## S = 993, p-value = 0.92
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.024781
cor.test(StrucMon$Leaf, StrucMon$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Leaf and StrucMon$Biomass
## S = 973, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0041301
cor.test(StrucMon$Leaf, StrucMon$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Leaf and StrucMon$Richness
## S = 908, p-value = 0.8
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.063115
cor.test(StrucMon$Leaf, StrucMon$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Leaf and StrucMon$Simpson_1.D
## S = 951, p-value = 0.94
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.01878
cor.test(StrucMon$Leaf, StrucMon$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Leaf and StrucMon$BasalHA
## S = 989, p-value = 0.94
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.020704

Branch Litter

cor.test(Struc$Branch, Struc$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$TreesHA
## S = 2907, p-value = 0.021
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.41399
cor.test(Struc$Branch, Struc$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$BasalHA
## S = 2981, p-value = 0.026
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.39903
cor.test(Struc$Branch, Struc$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$Biomass
## S = 3025, p-value = 0.03
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.39016
cor.test(Struc$Branch, Struc$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$Richness
## S = 3359, p-value = 0.077
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.32274
cor.test(Struc$Branch, Struc$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$Simpson_1.D
## S = 3927, p-value = 0.26
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.2083

Next we will use ANCOVAs to determine if landscapes differ in their relationship between Branch litter and the structural variables.

anova(lm(sqrt(Branch) ~ TreesHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                   Df Sum Sq Mean Sq F value Pr(>F)   
## TreesHA            1  0.165  0.1649    8.20  0.008 **
## Landscape          1  0.006  0.0065    0.32  0.574   
## TreesHA:Landscape  1  0.067  0.0670    3.33  0.079 . 
## Residuals         27  0.543  0.0201                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ BasalHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## BasalHA            1  0.117  0.1170    5.47  0.027 *
## Landscape          1  0.027  0.0268    1.26  0.272  
## BasalHA:Landscape  1  0.060  0.0601    2.81  0.105  
## Residuals         27  0.577  0.0214                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ Biomass * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## Biomass            1  0.109  0.1085    4.97  0.034 *
## Landscape          1  0.031  0.0310    1.42  0.244  
## Biomass:Landscape  1  0.052  0.0518    2.37  0.135  
## Residuals         27  0.590  0.0218                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ Richness * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                    Df Sum Sq Mean Sq F value Pr(>F)  
## Richness            1  0.093  0.0925    4.22   0.05 *
## Landscape           1  0.036  0.0356    1.63   0.21  
## Richness:Landscape  1  0.061  0.0608    2.77   0.11  
## Residuals          27  0.592  0.0219                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ Simpson_1.D * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                       Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D            1  0.048  0.0484    1.97   0.17
## Landscape              1  0.053  0.0531    2.16   0.15
## Simpson_1.D:Landscape  1  0.015  0.0146    0.59   0.45
## Residuals             27  0.665  0.0246
cor.test(StrucLome$Branch, StrucLome$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$TreesHA
## S = 72, p-value = 0.0016
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.8022
cor.test(StrucLome$Branch, StrucLome$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$Biomass
## S = 56, p-value = 0.00044
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.84615
cor.test(StrucLome$Branch, StrucLome$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$Richness
## S = 77.6, p-value = 0.0014
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.7868
cor.test(StrucLome$Branch, StrucLome$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$Simpson_1.D
## S = 140, p-value = 0.029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.61538
cor.test(StrucLome$Branch, StrucLome$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$BasalHA
## S = 60, p-value = 0.00065
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.83516
cor.test(StrucMon$Branch, StrucMon$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$TreesHA
## S = 861, p-value = 0.66
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.11151
cor.test(StrucMon$Branch, StrucMon$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$Biomass
## S = 793, p-value = 0.47
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.18172
cor.test(StrucMon$Branch, StrucMon$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$Richness
## S = 898, p-value = 0.77
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.073461
cor.test(StrucMon$Branch, StrucMon$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$Simpson_1.D
## S = 916, p-value = 0.83
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.054255
cor.test(StrucMon$Branch, StrucMon$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$BasalHA
## S = 830, p-value = 0.57
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.14389

Detritus Litter

cor.test(Struc$Detritus, Struc$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Detritus and Struc$TreesHA
## S = 3604, p-value = 0.14
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.27344
cor.test(Struc$Detritus, Struc$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Detritus and Struc$BasalHA
## S = 3014, p-value = 0.029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.39237
cor.test(Struc$Detritus, Struc$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Detritus and Struc$Biomass
## S = 3035, p-value = 0.031
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.38814
cor.test(Struc$Detritus, Struc$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Detritus and Struc$Richness
## S = 3369, p-value = 0.079
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.32072
cor.test(Struc$Detritus, Struc$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Detritus and Struc$Simpson_1.D
## S = 3417, p-value = 0.088
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.31114

Next we will use ANCOVAs to determine if landscapes differ in their relationship between Detritus litter and the structural variables.

anova(lm(sqrt(Detritus) ~ TreesHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## TreesHA            1  0.169  0.1690    3.57  0.070 .
## Landscape          1  0.020  0.0197    0.41  0.525  
## TreesHA:Landscape  1  0.172  0.1716    3.62  0.068 .
## Residuals         27  1.280  0.0474                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Detritus) ~ BasalHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## BasalHA            1  0.156  0.1556    3.14  0.088 .
## Landscape          1  0.003  0.0027    0.05  0.816  
## BasalHA:Landscape  1  0.145  0.1446    2.92  0.099 .
## Residuals         27  1.337  0.0495                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Detritus) ~ Biomass * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## Biomass            1  0.152  0.1519    3.02  0.094 .
## Landscape          1  0.001  0.0015    0.03  0.866  
## Biomass:Landscape  1  0.128  0.1275    2.53  0.123  
## Residuals         27  1.359  0.0503                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Detritus) ~ Richness * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus)
##                    Df Sum Sq Mean Sq F value Pr(>F)
## Richness            1  0.132  0.1321    2.57   0.12
## Landscape           1  0.001  0.0005    0.01   0.92
## Richness:Landscape  1  0.121  0.1213    2.36   0.14
## Residuals          27  1.386  0.0513
anova(lm(sqrt(Detritus) ~ Simpson_1.D * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Detritus)
##                       Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D            1  0.132  0.1318    2.39   0.13
## Landscape              1  0.001  0.0014    0.03   0.88
## Simpson_1.D:Landscape  1  0.017  0.0166    0.30   0.59
## Residuals             27  1.490  0.0552
cor.test(StrucLome$Detritus, StrucLome$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Detritus and StrucLome$TreesHA
## S = 82, p-value = 0.0029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.77473
cor.test(StrucLome$Detritus, StrucLome$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Detritus and StrucLome$Biomass
## S = 82, p-value = 0.0029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.77473
cor.test(StrucLome$Detritus, StrucLome$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Detritus and StrucLome$Richness
## S = 78.6, p-value = 0.0015
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.78404
cor.test(StrucLome$Detritus, StrucLome$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Detritus and StrucLome$Simpson_1.D
## S = 96, p-value = 0.0058
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.73626
cor.test(StrucLome$Detritus, StrucLome$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Detritus and StrucLome$BasalHA
## S = 76, p-value = 0.0021
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.79121
cor.test(StrucMon$Detritus, StrucMon$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Detritus and StrucMon$TreesHA
## S = 1053, p-value = 0.73
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.086732
cor.test(StrucMon$Detritus, StrucMon$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Detritus and StrucMon$Biomass
## S = 847, p-value = 0.62
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.12597
cor.test(StrucMon$Detritus, StrucMon$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Detritus and StrucMon$Richness
## S = 971, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0020693
cor.test(StrucMon$Detritus, StrucMon$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Detritus and StrucMon$Simpson_1.D
## S = 1007, p-value = 0.88
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.039648
cor.test(StrucMon$Detritus, StrucMon$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Detritus and StrucMon$BasalHA
## S = 885, p-value = 0.73
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.086957

Branch Litter

cor.test(Struc$Branch, Struc$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$TreesHA
## S = 2907, p-value = 0.021
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.41399
cor.test(Struc$Branch, Struc$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$BasalHA
## S = 2981, p-value = 0.026
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.39903
cor.test(Struc$Branch, Struc$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$Biomass
## S = 3025, p-value = 0.03
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.39016
cor.test(Struc$Branch, Struc$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$Richness
## S = 3359, p-value = 0.077
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.32274
cor.test(Struc$Branch, Struc$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Branch and Struc$Simpson_1.D
## S = 3927, p-value = 0.26
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.2083

Next we will use ANCOVAs to determine if landscapes differ in their relationship between Branch litter and the structural variables.

anova(lm(sqrt(Branch) ~ TreesHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                   Df Sum Sq Mean Sq F value Pr(>F)   
## TreesHA            1  0.165  0.1649    8.20  0.008 **
## Landscape          1  0.006  0.0065    0.32  0.574   
## TreesHA:Landscape  1  0.067  0.0670    3.33  0.079 . 
## Residuals         27  0.543  0.0201                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ BasalHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## BasalHA            1  0.117  0.1170    5.47  0.027 *
## Landscape          1  0.027  0.0268    1.26  0.272  
## BasalHA:Landscape  1  0.060  0.0601    2.81  0.105  
## Residuals         27  0.577  0.0214                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ Biomass * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## Biomass            1  0.109  0.1085    4.97  0.034 *
## Landscape          1  0.031  0.0310    1.42  0.244  
## Biomass:Landscape  1  0.052  0.0518    2.37  0.135  
## Residuals         27  0.590  0.0218                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ Richness * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                    Df Sum Sq Mean Sq F value Pr(>F)  
## Richness            1  0.093  0.0925    4.22   0.05 *
## Landscape           1  0.036  0.0356    1.63   0.21  
## Richness:Landscape  1  0.061  0.0608    2.77   0.11  
## Residuals          27  0.592  0.0219                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Branch) ~ Simpson_1.D * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Branch)
##                       Df Sum Sq Mean Sq F value Pr(>F)
## Simpson_1.D            1  0.048  0.0484    1.97   0.17
## Landscape              1  0.053  0.0531    2.16   0.15
## Simpson_1.D:Landscape  1  0.015  0.0146    0.59   0.45
## Residuals             27  0.665  0.0246
cor.test(StrucLome$Branch, StrucLome$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$TreesHA
## S = 72, p-value = 0.0016
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.8022
cor.test(StrucLome$Branch, StrucLome$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$Biomass
## S = 56, p-value = 0.00044
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.84615
cor.test(StrucLome$Branch, StrucLome$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$Richness
## S = 77.6, p-value = 0.0014
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.7868
cor.test(StrucLome$Branch, StrucLome$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$Simpson_1.D
## S = 140, p-value = 0.029
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.61538
cor.test(StrucLome$Branch, StrucLome$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Branch and StrucLome$BasalHA
## S = 60, p-value = 0.00065
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.83516
cor.test(StrucMon$Branch, StrucMon$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$TreesHA
## S = 861, p-value = 0.66
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.11151
cor.test(StrucMon$Branch, StrucMon$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$Biomass
## S = 793, p-value = 0.47
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.18172
cor.test(StrucMon$Branch, StrucMon$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$Richness
## S = 898, p-value = 0.77
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.073461
cor.test(StrucMon$Branch, StrucMon$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$Simpson_1.D
## S = 916, p-value = 0.83
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.054255
cor.test(StrucMon$Branch, StrucMon$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Branch and StrucMon$BasalHA
## S = 830, p-value = 0.57
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.14389

Flower Litter

cor.test(Struc$Flower, Struc$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Flower and Struc$TreesHA
## S = 3216, p-value = 0.052
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.35153
cor.test(Struc$Flower, Struc$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Flower and Struc$BasalHA
## S = 3102, p-value = 0.038
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.37468
cor.test(Struc$Flower, Struc$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Flower and Struc$Biomass
## S = 3219, p-value = 0.053
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.35109
cor.test(Struc$Flower, Struc$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Flower and Struc$Richness
## S = 2553, p-value = 0.0057
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.48526
cor.test(Struc$Flower, Struc$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Flower and Struc$Simpson_1.D
## S = 2573, p-value = 0.0061
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.4812

Next we will use ANCOVAs to determine if landscapes differ in their relationship between Flower litter and the structural variables.

anova(lm(sqrt(Flower) ~ TreesHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Flower)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## TreesHA            1 0.0236 0.02359    3.92  0.058 .
## Landscape          1 0.0090 0.00904    1.50  0.231  
## TreesHA:Landscape  1 0.0284 0.02840    4.72  0.039 *
## Residuals         27 0.1625 0.00602                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Flower) ~ BasalHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Flower)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## BasalHA            1 0.0220 0.02201    3.70  0.065 .
## Landscape          1 0.0147 0.01474    2.48  0.127  
## BasalHA:Landscape  1 0.0260 0.02602    4.37  0.046 *
## Residuals         27 0.1608 0.00595                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Flower) ~ Biomass * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Flower)
##                   Df Sum Sq Mean Sq F value Pr(>F)  
## Biomass            1 0.0207 0.02068    3.44  0.075 .
## Landscape          1 0.0160 0.01598    2.66  0.115  
## Biomass:Landscape  1 0.0246 0.02455    4.08  0.053 .
## Residuals         27 0.1623 0.00601                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Flower) ~ Richness * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Flower)
##                    Df Sum Sq Mean Sq F value Pr(>F)   
## Richness            1 0.0434  0.0434    7.96 0.0089 **
## Landscape           1 0.0164  0.0164    3.00 0.0944 . 
## Richness:Landscape  1 0.0166  0.0166    3.04 0.0928 . 
## Residuals          27 0.1472  0.0055                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Flower) ~ Simpson_1.D * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Flower)
##                       Df Sum Sq Mean Sq F value Pr(>F)  
## Simpson_1.D            1 0.0238 0.02383    3.70  0.065 .
## Landscape              1 0.0246 0.02461    3.82  0.061 .
## Simpson_1.D:Landscape  1 0.0012 0.00121    0.19  0.669  
## Residuals             27 0.1739 0.00644                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.test(StrucLome$Flower, StrucLome$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Flower and StrucLome$TreesHA
## S = 81.2, p-value = 0.0018
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.77686
cor.test(StrucLome$Flower, StrucLome$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Flower and StrucLome$Biomass
## S = 73.2, p-value = 0.0011
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.7989
cor.test(StrucLome$Flower, StrucLome$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Flower and StrucLome$Richness
## S = 58.7, p-value = 0.00034
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.83862
cor.test(StrucLome$Flower, StrucLome$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Flower and StrucLome$Simpson_1.D
## S = 91.2, p-value = 0.0032
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.74931
cor.test(StrucLome$Flower, StrucLome$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Flower and StrucLome$BasalHA
## S = 81.2, p-value = 0.0018
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.77686
cor.test(StrucMon$Flower, StrucMon$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Flower and StrucMon$TreesHA
## S = 977, p-value = 0.97
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0082687
cor.test(StrucMon$Flower, StrucMon$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Flower and StrucMon$Biomass
## S = 969, p-value = 1
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho 
##   0
cor.test(StrucMon$Flower, StrucMon$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Flower and StrucMon$Richness
## S = 695, p-value = 0.26
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.28276
cor.test(StrucMon$Flower, StrucMon$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Flower and StrucMon$Simpson_1.D
## S = 587, p-value = 0.11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.39375
cor.test(StrucMon$Flower, StrucMon$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Flower and StrucMon$BasalHA
## S = 952, p-value = 0.94
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.017617

Fruit Litter

cor.test(Struc$Fruit, Struc$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Fruit and Struc$TreesHA
## S = 5340, p-value = 0.68
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.076636
cor.test(Struc$Fruit, Struc$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Fruit and Struc$BasalHA
## S = 4587, p-value = 0.69
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.075255
cor.test(Struc$Fruit, Struc$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Fruit and Struc$Biomass
## S = 4495, p-value = 0.62
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.093769
cor.test(Struc$Fruit, Struc$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Fruit and Struc$Richness
## S = 4608, p-value = 0.7
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.071054
cor.test(Struc$Fruit, Struc$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Fruit and Struc$Simpson_1.D
## S = 4586, p-value = 0.69
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.075369

Next we will use ANCOVAs to determine if landscapes differ in their relationship between Fruit litter and the structural variables.

anova(lm(sqrt(Fruit) ~ TreesHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Fruit)
##                   Df Sum Sq Mean Sq F value Pr(>F)
## TreesHA            1  0.000 0.00027    0.01   0.91
## Landscape          1  0.029 0.02894    1.30   0.26
## TreesHA:Landscape  1  0.001 0.00089    0.04   0.84
## Residuals         27  0.602 0.02230
anova(lm(sqrt(Fruit) ~ BasalHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Fruit)
##                   Df Sum Sq Mean Sq F value Pr(>F)
## BasalHA            1  0.006 0.00554    0.25   0.62
## Landscape          1  0.028 0.02774    1.25   0.27
## BasalHA:Landscape  1  0.001 0.00143    0.06   0.80
## Residuals         27  0.598 0.02213
anova(lm(sqrt(Fruit) ~ Biomass * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Fruit)
##                   Df Sum Sq Mean Sq F value Pr(>F)
## Biomass            1  0.007 0.00707    0.32   0.58
## Landscape          1  0.027 0.02699    1.22   0.28
## Biomass:Landscape  1  0.000 0.00040    0.02   0.89
## Residuals         27  0.598 0.02214
anova(lm(sqrt(Fruit) ~ Richness * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Fruit)
##                    Df Sum Sq Mean Sq F value Pr(>F)
## Richness            1  0.003 0.00339    0.15   0.70
## Landscape           1  0.025 0.02543    1.14   0.29
## Richness:Landscape  1  0.001 0.00139    0.06   0.80
## Residuals          27  0.602 0.02230
anova(lm(sqrt(Fruit) ~ Simpson_1.D * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Fruit)
##                       Df Sum Sq Mean Sq F value Pr(>F)  
## Simpson_1.D            1  0.080  0.0803    4.12  0.052 .
## Landscape              1  0.016  0.0157    0.80  0.378  
## Simpson_1.D:Landscape  1  0.010  0.0098    0.50  0.485  
## Residuals             27  0.526  0.0195                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.test(StrucLome$Fruit, StrucLome$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Fruit and StrucLome$TreesHA
## S = 366, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0054945
cor.test(StrucLome$Fruit, StrucLome$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Fruit and StrucLome$Biomass
## S = 332, p-value = 0.78
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.087912
cor.test(StrucLome$Fruit, StrucLome$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Fruit and StrucLome$Richness
## S = 306, p-value = 0.6
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.15956
cor.test(StrucLome$Fruit, StrucLome$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Fruit and StrucLome$Simpson_1.D
## S = 348, p-value = 0.89
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.043956
cor.test(StrucLome$Fruit, StrucLome$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Fruit and StrucLome$BasalHA
## S = 334, p-value = 0.79
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.082418
cor.test(StrucMon$Fruit, StrucMon$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Fruit and StrucMon$TreesHA
## S = 1098, p-value = 0.6
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## -0.1332
cor.test(StrucMon$Fruit, StrucMon$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Fruit and StrucMon$Biomass
## S = 911, p-value = 0.81
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.059886
cor.test(StrucMon$Fruit, StrucMon$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Fruit and StrucMon$Richness
## S = 982, p-value = 0.96
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.013451
cor.test(StrucMon$Fruit, StrucMon$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Fruit and StrucMon$Simpson_1.D
## S = 976, p-value = 0.98
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0073035
cor.test(StrucMon$Fruit, StrucMon$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Fruit and StrucMon$BasalHA
## S = 946, p-value = 0.93
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.02381

Total Litter

cor.test(Struc$Total, Struc$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Total and Struc$TreesHA
## S = 4572, p-value = 0.68
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.078242
cor.test(Struc$Total, Struc$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Total and Struc$BasalHA
## S = 3535, p-value = 0.12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.28727
cor.test(Struc$Total, Struc$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Total and Struc$Biomass
## S = 3501, p-value = 0.11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.29418
cor.test(Struc$Total, Struc$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Total and Struc$Richness
## S = 3585, p-value = 0.13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.27712
cor.test(Struc$Total, Struc$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  Struc$Total and Struc$Simpson_1.D
## S = 3392, p-value = 0.083
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.31619

Next we will use ANCOVAs to determine if landscapes differ in their relationship between Total litter and the structural variables.

anova(lm(sqrt(Total) ~ TreesHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Total)
##                   Df Sum Sq Mean Sq F value  Pr(>F)    
## TreesHA            1  0.275   0.275    2.99 0.09513 .  
## Landscape          1  1.375   1.375   14.97 0.00062 ***
## TreesHA:Landscape  1  0.816   0.816    8.89 0.00602 ** 
## Residuals         27  2.480   0.092                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Total) ~ BasalHA * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Total)
##                   Df Sum Sq Mean Sq F value Pr(>F)   
## BasalHA            1  0.442   0.442    4.35 0.0466 * 
## Landscape          1  1.075   1.075   10.58 0.0031 **
## BasalHA:Landscape  1  0.685   0.685    6.74 0.0150 * 
## Residuals         27  2.744   0.102                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Total) ~ Biomass * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Total)
##                   Df Sum Sq Mean Sq F value Pr(>F)   
## Biomass            1  0.486   0.486    4.63 0.0405 * 
## Landscape          1  1.027   1.027    9.79 0.0042 **
## Biomass:Landscape  1  0.602   0.602    5.74 0.0237 * 
## Residuals         27  2.831   0.105                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Total) ~ Richness * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Total)
##                    Df Sum Sq Mean Sq F value Pr(>F)   
## Richness            1  0.507   0.507    4.67 0.0396 * 
## Landscape           1  0.973   0.973    8.97 0.0058 **
## Richness:Landscape  1  0.537   0.537    4.95 0.0346 * 
## Residuals          27  2.929   0.108                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(sqrt(Total) ~ Simpson_1.D * Landscape, data = Struc)) 
## Analysis of Variance Table
## 
## Response: sqrt(Total)
##                       Df Sum Sq Mean Sq F value Pr(>F)  
## Simpson_1.D            1   0.44   0.441    3.23  0.083 .
## Landscape              1   0.77   0.773    5.66  0.025 *
## Simpson_1.D:Landscape  1   0.04   0.045    0.33  0.573  
## Residuals             27   3.69   0.137                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cor.test(StrucLome$Total, StrucLome$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Total and StrucLome$TreesHA
## S = 56, p-value = 0.00044
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.84615
cor.test(StrucLome$Total, StrucLome$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Total and StrucLome$Biomass
## S = 58, p-value = 0.00054
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.84066
cor.test(StrucLome$Total, StrucLome$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Total and StrucLome$Richness
## S = 52.6, p-value = 0.00019
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.85557
cor.test(StrucLome$Total, StrucLome$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Total and StrucLome$Simpson_1.D
## S = 86, p-value = 0.0036
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.76374
cor.test(StrucLome$Total, StrucLome$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucLome$Total and StrucLome$BasalHA
## S = 50, p-value = 0.00019
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.86264
cor.test(StrucMon$Total, StrucMon$TreesHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Total and StrucMon$TreesHA
## S = 1120, p-value = 0.54
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.15591
cor.test(StrucMon$Total, StrucMon$Biomass, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Total and StrucMon$Biomass
## S = 972, p-value = 0.99
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0030976
cor.test(StrucMon$Total, StrucMon$Richness, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Total and StrucMon$Richness
## S = 1041, p-value = 0.77
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.074496
cor.test(StrucMon$Total, StrucMon$Simpson_1.D, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Total and StrucMon$Simpson_1.D
## S = 1092, p-value = 0.61
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## -0.12729
cor.test(StrucMon$Total, StrucMon$BasalHA, method="spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  StrucMon$Total and StrucMon$BasalHA
## S = 1020, p-value = 0.84
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.052795

Figure 5

Let’s use Ordination techniques to evaluate how plots differentiate. First we make a copy of Struc and call it StrucNew StrucNew is then edited to rename the factor levels of the Succession variable. We then subset the columns we don’t need, rename StrucNew as PCAStrucNew, use prcomp to draft the ordination, rename our landscapes using the list function, relevel using the relevel function, and graph using autoplot.

StrucNew <- Struc
levels(StrucNew$Succession) <- list(Early  = "DM1", Intermediate = "DM2", `Old-Growth` = "RM", 
                                    Early  = "DL1", Intermediate = "DL2", `Old-Growth` = "RL")
StrucNew2 <- StrucNew[5:19]
StrucNew3 <- StrucNew2[ , -which(names(StrucNew2) %in% c("Shannon_H","Chao.1", "Family", "Biomass", "Total.Mg.ha.month", "Leaf", "Branch", "Flower", "Fruit", "Detritus"))]

PCAStrucNew <- StrucNew

pca_resx <- prcomp(StrucNew3, scale. = TRUE)
levels(PCAStrucNew$Landscape) <- list(Mountain  = "Mountain", Hill = "Lomerio") 
PCAStrucNew$Landscape <- relevel(PCAStrucNew$Landscape, ref = "Hill") 

autoplot(pca_resx, data = PCAStrucNew, colour = "Landscape", shape = "Succession",
                    loadings = TRUE, loadings.colour = c('#AA54B4', '#AA54B4', '#AA54B4',
                                                         '#AA54B4', 
                                                         '#B0B315', '#3E99D9'),
                    loadings.label = TRUE, loadings.label.size = 4, frame = TRUE, frame.type = "norm", 
                    loadings.label.label = c('Tree Density', 'Basal Area', 'Species Density',
                                             'Diversity', 
                                             'Total Litter',  'Precipitation'),
                    loadings.label.colour = c('#AA54B4', '#AA54B4', '#AA54B4',
                                              '#AA54B4', 
                                              '#B0B315', '#3E99D9'),
                    loadings.label.repel=T, loadings.label.hjust = -0.35,
                    loadings.label.vjust = 0.25) +
  guides(fill = "none") +
  scale_colour_manual(values=c("darkgreen", "brown")) + 
  scale_fill_manual(values=c("#23431F", "#67360D")) + 
  theme_bw(base_size = 15) + guides(color=guide_legend(title="")) + 
  guides(shape=guide_legend(title="")) +
  theme(legend.position = "bottom", 
        legend.box = "vertical",
        legend.spacing.y = unit(0, "pt"), 
        legend.margin = margin(t = 1, b = 1))

Next, we use the cor function within the eigen function to obtain our Eigvenvalues. We then do some more magic to get in the format we need!

Eigenvalues

eigen(cor(StrucNew3))
## eigen() decomposition
## $values
## [1] 3.457685 1.527042 0.520299 0.337626 0.110986 0.046362
## 
## $vectors
##           [,1]       [,2]      [,3]      [,4]      [,5]       [,6]
## [1,] -0.497313  0.1774644  0.229850  0.080832  0.750232  0.3145994
## [2,] -0.518787  0.0588992  0.137140  0.267637 -0.096461 -0.7922433
## [3,] -0.514155  0.0058209  0.044145  0.304629 -0.607277  0.5216107
## [4,] -0.408332 -0.1444513 -0.805223 -0.399123  0.063821 -0.0253393
## [5,] -0.208078 -0.6386053  0.488023 -0.552321 -0.075126 -0.0041811
## [6,]  0.095666 -0.7323390 -0.199654  0.603888  0.222138  0.0253084
EigenTable2A <- as.data.frame(eigen(cor(StrucNew3))$values)
EigenTable2A <- as.data.frame(EigenTable2A)
colnames(EigenTable2A) <- "Eigenvalue"

EigenTable2A$CumulativePercentVariation <- EigenTable2A$Eigenvalue/sum(EigenTable2A$Eigenvalue)*100
EigenTable2A <- t(EigenTable2A) # "t" function transposes. 

EigenTable2B <- as.data.frame(eigen(cor(StrucNew3))$vectors)
EigenTable <- rbind(EigenTable2A,EigenTable2B)

row.names(EigenTable) <-c("Eigenvalue", "Cumulative Percent Variation", colnames(StrucNew3))
colnames(EigenTable) <- c("PC1", "PC2", 
                          "PC3", "PC4", "PC5", "PC6")

kbl(EigenTable) %>% kable_classic()
PC1 PC2 PC3 PC4 PC5 PC6
Eigenvalue 3.45769 1.52704 0.52030 0.33763 0.11099 0.04636
Cumulative Percent Variation 57.62809 25.45070 8.67165 5.62709 1.84977 0.77270
TreesHA -0.49731 0.17746 0.22985 0.08083 0.75023 0.31460
BasalHA -0.51879 0.05890 0.13714 0.26764 -0.09646 -0.79224
Richness -0.51416 0.00582 0.04415 0.30463 -0.60728 0.52161
Simpson_1.D -0.40833 -0.14445 -0.80522 -0.39912 0.06382 -0.02534
Total -0.20808 -0.63861 0.48802 -0.55232 -0.07513 -0.00418
RainNow 0.09567 -0.73234 -0.19965 0.60389 0.22214 0.02531

Structure x Habitat

Do habitats within landscapes differ in structure? We can find out via anova function on an lm function of each variable by Age Group.

anova(lm(TreesHA ~ Succession, data = StrucLome))
## Analysis of Variance Table
## 
## Response: TreesHA
##            Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession  2 344476  172238    48.3 7.3e-06 ***
## Residuals  10  35669    3567                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(BasalHA ~ Succession, data = StrucLome))
## Analysis of Variance Table
## 
## Response: BasalHA
##            Df Sum Sq Mean Sq F value Pr(>F)    
## Succession  2    865     432    58.6  3e-06 ***
## Residuals  10     74       7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Simpson_1.D ~ Succession, data = StrucLome))
## Analysis of Variance Table
## 
## Response: Simpson_1.D
##            Df Sum Sq Mean Sq F value Pr(>F)  
## Succession  2  0.139  0.0694    3.19  0.085 .
## Residuals  10  0.218  0.0218                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(Richness ~ Succession, data = StrucLome))
## Analysis of Variance Table
## 
## Response: Richness
##            Df Sum Sq Mean Sq F value  Pr(>F)    
## Succession  2   8747    4374    62.4 2.2e-06 ***
## Residuals  10    701      70                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lm(TreesHA ~ Succession, data = StrucMon))
## Analysis of Variance Table
## 
## Response: TreesHA
##            Df  Sum Sq Mean Sq F value Pr(>F)
## Succession  2   11916    5958    0.08   0.92
## Residuals  15 1071344   71423
anova(lm(BasalHA ~ Succession, data = StrucMon))
## Analysis of Variance Table
## 
## Response: BasalHA
##            Df Sum Sq Mean Sq F value Pr(>F)
## Succession  2    125    62.4    0.52    0.6
## Residuals  15   1789   119.3
anova(lm(Simpson_1.D ~ Succession, data = StrucMon))
## Analysis of Variance Table
## 
## Response: Simpson_1.D
##            Df Sum Sq Mean Sq F value Pr(>F)
## Succession  2  0.060  0.0302    0.93   0.42
## Residuals  15  0.488  0.0325
anova(lm(Richness ~ Succession, data = StrucMon))
## Analysis of Variance Table
## 
## Response: Richness
##            Df Sum Sq Mean Sq F value Pr(>F)
## Succession  2   1272     636    0.58   0.57
## Residuals  15  16591    1106

CONCLUSION: Differences in Hill, not in Mountain! Let’s get some descriptive statistics for the structural variables now.

kbl(summarySE(StrucLome, measurevar="TreesHA", groupvars=c("Succession"), 
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N TreesHA se
DL1 3 148.00 45.078
DL2 6 446.67 26.385
RL 4 592.00 14.697
kbl(summarySE(StrucLome, measurevar="BasalHA", groupvars=c("Succession"), 
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N BasalHA se
DL1 3 3.3733 1.66455
DL2 6 13.2167 1.35706
RL 4 25.5000 0.39791
kbl(summarySE(StrucLome, measurevar="Simpson_1.D", groupvars=c("Succession"),
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N Simpson_1.D se
DL1 3 0.70800 0.18804
DL2 6 0.93633 0.01339
RL 4 0.97200 0.00280
kbl(summarySE(StrucLome, measurevar="Richness", groupvars=c("Succession"),
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N Richness se
DL1 3 10.667 3.1798
DL2 6 38.667 2.9059
RL 4 80.250 5.6771
kbl(summarySE(StrucMon, measurevar="TreesHA", groupvars=c("Succession"), 
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N TreesHA se
DM1 6 547.33 132.111
DM2 6 553.33 77.015
RM 6 604.67 111.026
kbl(summarySE(StrucMon, measurevar="BasalHA", groupvars=c("Succession"), 
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N BasalHA se
DM1 6 15.167 4.5491
DM2 6 15.333 3.7830
RM 6 20.833 4.9626
kbl(summarySE(StrucMon, measurevar="Simpson_1.D", groupvars=c("Succession"), 
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N Simpson_1.D se
DM1 6 0.82500 0.09283
DM2 6 0.80500 0.08586
RM 6 0.93667 0.01687
kbl(summarySE(StrucMon, measurevar="Richness", groupvars=c("Succession"), 
    na.rm = TRUE)[,-4][,-5][,-7], format = "html", 
    table.attr = "style='width:50%;'", row.names = FALSE) %>% 
  kable_classic()
Succession N Richness se
DM1 6 42.333 13.696
DM2 6 42.000 13.760
RM 6 60.000 13.272

Species Composition

We can look at which were the dominant species by examineing the Composition dataframe. First, we will need to join DP and EF as that’s the classification in the MS. Next, we subset so that we can obtain Hill and Mountain dataframes.

levels(Composition$Chronosecuence) <- list("EF" = "DP",  "EF" = "EF", "IF"  = "IF", "OF" = "OF") 

SPHill <- droplevels(subset(Composition, Composition$Landscape == "Hill"))
SPMountain <- droplevels(subset(Composition, Composition$Landscape == "Mountain"))

Next we can use the aggregate function to sum up the abundance by Species and Chronosequence. Next we order the dataframe by abundance.

SPCountH <- as.data.frame(aggregate(abun ~ Species + Chronosecuence, data = SPHill, FUN = sum, na.rm = TRUE)) 
SPCountH <-SPCountH[order(SPCountH$abun, decreasing = TRUE),]
SPCountH[1:15,]
##                         Species Chronosecuence abun
## 172         Siparuna guianensis             EF  242
## 111               Miconia elata             EF  160
## 441         Siparuna guianensis             IF  119
## 197   Adenocalymma cladotrichum             IF  114
## 459         Tapirira guianensis             IF  104
## 98               Isertia laevis             EF   76
## 2     Adenocalymma cladotrichum             EF   74
## 421 Pseudosenefeldera inclinata             IF   71
## 80      Henriettea fascicularis             EF   69
## 149          Piptocoma discolor             EF   65
## 687 Pseudosenefeldera inclinata             OF   65
## 306     Henriettea fascicularis             IF   59
## 112         Miconia minutiflora             EF   55
## 620         Miconia pterocaulon             OF   53
## 294          Guatteria punctata             IF   49

Great, now we can subset by Age Group within hill and take a look.

SPCountH.EF <- droplevels(subset(SPCountH, SPCountH$Chronosecuence == "EF"))
SPCountH.IF <- droplevels(subset(SPCountH, SPCountH$Chronosecuence == "IF"))
SPCountH.OF <- droplevels(subset(SPCountH, SPCountH$Chronosecuence == "OF"))
SPCountH.EF[1:5,]; SPCountH.IF[1:5,]; SPCountH.OF[1:5,]
##                       Species Chronosecuence abun
## 172       Siparuna guianensis             EF  242
## 111             Miconia elata             EF  160
## 98             Isertia laevis             EF   76
## 2   Adenocalymma cladotrichum             EF   74
## 80    Henriettea fascicularis             EF   69
##                         Species Chronosecuence abun
## 441         Siparuna guianensis             IF  119
## 197   Adenocalymma cladotrichum             IF  114
## 459         Tapirira guianensis             IF  104
## 421 Pseudosenefeldera inclinata             IF   71
## 306     Henriettea fascicularis             IF   59
##                         Species Chronosecuence abun
## 687 Pseudosenefeldera inclinata             OF   65
## 620         Miconia pterocaulon             OF   53
## 729             Virola elongata             OF   40
## 623           Miconia splendens             OF   37
## 649           Oenocarpus bataua             OF   27

And now we do the same for Mountain.

SPCountM <- as.data.frame(aggregate(abun ~ Species + Chronosecuence, 
                                    data = SPMountain, FUN = sum, na.rm = TRUE)) 
SPCountM <-SPCountM[order(SPCountM$abun, decreasing = TRUE),]
SPCountM[1:15,]
##                          Species Chronosecuence abun
## 64       Henriettea fascicularis             EF  459
## 936  Pseudosenefeldera inclinata             OF  168
## 336  Graffenrieda conostegioides             IF  163
## 1020          Wettinia praemorsa             OF  130
## 106          Miconia minutiflora             EF  114
## 104          Miconia lourtegiana             EF  103
## 1014             Virola elongata             OF   73
## 47           Eugenia lambertiana             EF   68
## 80              Inga thibaudiana             EF   66
## 187             Vismia baccifera             EF   66
## 240             Casearia arborea             IF   62
## 121              Myrcia splenden             EF   60
## 795       Ladenbergia muzonensis             OF   56
## 16      Bellucia grossularioides             EF   52
## 736      Graffenrieda colombiana             OF   50
SPCountM.EF <- droplevels(subset(SPCountM, SPCountM$Chronosecuence == "EF"))
SPCountM.IF <- droplevels(subset(SPCountM, SPCountM$Chronosecuence == "IF"))
SPCountM.OF <- droplevels(subset(SPCountM, SPCountM$Chronosecuence == "OF"))
SPCountM.EF[1:5,]; SPCountM.IF[1:5,]; SPCountM.OF[1:5,]
##                     Species Chronosecuence abun
## 64  Henriettea fascicularis             EF  459
## 106     Miconia minutiflora             EF  114
## 104     Miconia lourtegiana             EF  103
## 47      Eugenia lambertiana             EF   68
## 80         Inga thibaudiana             EF   66
##                         Species Chronosecuence abun
## 336 Graffenrieda conostegioides             IF  163
## 240            Casearia arborea             IF   62
## 418           Matayba inelegans             IF   46
## 279           Cyathea lasiosora             IF   45
## 394    Ladenbergia oblongifolia             IF   41
##                          Species Chronosecuence abun
## 936  Pseudosenefeldera inclinata             OF  168
## 1020          Wettinia praemorsa             OF  130
## 1014             Virola elongata             OF   73
## 795       Ladenbergia muzonensis             OF   56
## 736      Graffenrieda colombiana             OF   50

Regressions

Let’s calculate regressions for Stand Age vs. Leaf Litter Production. First we need to manually insert the ages of the plots. These must be inserted in this order.

Struc$Age <- as.numeric(c("12", "7", "10", "25", "35", "30", "28", "30", "30", 
            "12", "5", "5", "17", "5", "13", "26", "35", "31", "37", "22", 
            "24", "40", "57", "49", "55", "47", "55", "60", "65", "68", "49"))

Next we relevel and rename.

levels(Struc$Landscape) <- list(Hill  = "Lomerio", Mountain = "Mountain")
Struc$`Stand Age Group` <- Struc$Succession
levels(Struc$`Stand Age Group`) <- list(EH  = "DL1", IH = "DL2", OH = "RL",
                                        EM = "DM1", IM = "DM2", OM = "RM")

Figure 3

We’re now ready for graphing! For this we will use ggplot and ggarrange.

GGNew1 <- ggplot(Struc, aes(x=Age, y=Leaf, group = Landscape, color = Landscape)) +
  geom_point(size=1, alpha=0.9) +
  geom_smooth(size = 0.5, method=lm, aes(fill = Landscape), color = "black") +
  scale_color_manual(values=c("#82C27B", "#FBBA84")) + 
  scale_fill_manual(values=c("#82C27B", "#FBBA84")) + 
  ylab(bquote('Leaf Litter Production '(Mg~ha^-1~mth^-1))) + 
  xlab("Stand Age") + 
  theme_bw(base_size = 9) +
  guides(color = "none") + 
  geom_vline(xintercept = c(20, 40), linetype = "dashed", color = "black")

GGNew2 <- ggplot(Struc, aes(x=Age, y=Total, group = Landscape , color = Landscape)) +
  geom_point(size=1, alpha=0.9) +
  geom_smooth(size = 0.5, method=lm, aes(fill = Landscape), color = "black") +
  scale_color_manual(values=c("#82C27B", "#FBBA84")) + 
  scale_fill_manual(values=c("#82C27B", "#FBBA84")) + 
  ylab(bquote('Total Litter Production '(Mg~ha^-1~mth^-1))) + 
  xlab("Stand Age") + 
  theme_bw(base_size = 9) +
  guides(color = "none") + 
  geom_vline(xintercept = c(20, 40), linetype = "dashed", color = "black")

ggarrange(GGNew2, GGNew1, 
          align='h', labels = "auto", legend = "right", 
          common.legend = T, nrow = 2, ncol = 1)

Session Report

## Analyses were conducted using the R Statistical language (version 4.2.2; R Core
## Team, 2022) on macOS Ventura 13.2.1, using the packages ggthemes (version
## 4.2.4; Arnold J, 2021), lme4 (version 1.1.31; Bates D et al., 2015), Matrix
## (version 1.5.1; Bates D et al., 2022), conover.test (version 1.1.5; Dinno A,
## 2017), lubridate (version 1.9.2; Grolemund G, Wickham H, 2011), Rmisc (version
## 1.5.1; Hope RM, 2022), ggpubr (version 0.6.0; Kassambara A, 2023), lmerTest
## (version 3.1.3; Kuznetsova A et al., 2017), emmeans (version 1.8.4.1; Lenth R,
## 2023), report (version 0.5.6; Makowski D et al., 2023), patchwork (version
## 1.1.2; Pedersen T, 2022), lattice (version 0.20.45; Sarkar D, 2008), ggfortify
## (version 0.4.15; Tang Y et al., 2016), reshape (version 0.8.9; Wickham H,
## 2007), plyr (version 1.8.8; Wickham H, 2011), ggplot2 (version 3.4.1; Wickham
## H, 2016), dplyr (version 1.1.0; Wickham H et al., 2023), zoo (version 1.8.11;
## Zeileis A, Grothendieck G, 2005), lmtest (version 0.9.40; Zeileis A, Hothorn T,
## 2002) and kableExtra (version 1.3.4; Zhu H, 2021).
## 
## References
## ----------
##   - Arnold J (2021). _ggthemes: Extra Themes, Scales and Geoms for 'ggplot2'_. R
## package version 4.2.4, <https://CRAN.R-project.org/package=ggthemes>.
##   - Bates D, Mächler M, Bolker B, Walker S (2015). "Fitting Linear Mixed-Effects
## Models Using lme4." _Journal of Statistical Software_, *67*(1), 1-48.
## doi:10.18637/jss.v067.i01 <https://doi.org/10.18637/jss.v067.i01>.
##   - Bates D, Maechler M, Jagan M (2022). _Matrix: Sparse and Dense Matrix Classes
## and Methods_. R package version 1.5-1,
## <https://CRAN.R-project.org/package=Matrix>.
##   - Dinno A (2017). _conover.test: Conover-Iman Test of Multiple Comparisons
## Using Rank Sums_. R package version 1.1.5,
## <https://CRAN.R-project.org/package=conover.test>.
##   - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate."
## _Journal of Statistical Software_, *40*(3), 1-25.
## <https://www.jstatsoft.org/v40/i03/>.
##   - Hope RM (2022). _Rmisc: Ryan Miscellaneous_. R package version 1.5.1,
## <https://CRAN.R-project.org/package=Rmisc>.
##   - Kassambara A (2023). _ggpubr: 'ggplot2' Based Publication Ready Plots_. R
## package version 0.6.0, <https://CRAN.R-project.org/package=ggpubr>.
##   - Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package: Tests
## in Linear Mixed Effects Models." _Journal of Statistical Software_, *82*(13),
## 1-26. doi:10.18637/jss.v082.i13 <https://doi.org/10.18637/jss.v082.i13>.
##   - Lenth R (2023). _emmeans: Estimated Marginal Means, aka Least-Squares Means_.
## R package version 1.8.4-1, <https://CRAN.R-project.org/package=emmeans>.
##   - Makowski D, Lüdecke D, Patil I, Thériault R (2023). "Automated Results
## Reporting as a Practical Tool to Improve Reproducibility and Methodological
## Best Practices Adoption." _CRAN_. <https://easystats.github.io/report/>.
##   - Pedersen T (2022). _patchwork: The Composer of Plots_. R package version
## 1.1.2, <https://CRAN.R-project.org/package=patchwork>.
##   - R Core Team (2022). _R: A Language and Environment for Statistical
## Computing_. R Foundation for Statistical Computing, Vienna, Austria.
## <https://www.R-project.org/>.
##   - Sarkar D (2008). _Lattice: Multivariate Data Visualization with R_. Springer,
## New York. ISBN 978-0-387-75968-5, <http://lmdvr.r-forge.r-project.org>.
##   - Tang Y, Horikoshi M, Li W (2016). "ggfortify: Unified Interface to Visualize
## Statistical Result of Popular R Packages." _The R Journal_, *8*(2), 474-485.
## doi:10.32614/RJ-2016-060 <https://doi.org/10.32614/RJ-2016-060>,
## <https://doi.org/10.32614/RJ-2016-060>. Horikoshi M, Tang Y (2018). _ggfortify:
## Data Visualization Tools for Statistical Analysis Results_.
## <https://CRAN.R-project.org/package=ggfortify>.
##   - Wickham H (2007). "Reshaping data with the reshape package." _Journal of
## Statistical Software_, *21*(12). <https://www.jstatsoft.org/v21/i12/>.
##   - Wickham H (2011). "The Split-Apply-Combine Strategy for Data Analysis."
## _Journal of Statistical Software_, *40*(1), 1-29.
## <https://www.jstatsoft.org/v40/i01/>.
##   - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_.
## Springer-Verlag New York. ISBN 978-3-319-24277-4,
## <https://ggplot2.tidyverse.org>.
##   - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar
## of Data Manipulation_. R package version 1.1.0,
## <https://CRAN.R-project.org/package=dplyr>.
##   - Zeileis A, Grothendieck G (2005). "zoo: S3 Infrastructure for Regular and
## Irregular Time Series." _Journal of Statistical Software_, *14*(6), 1-27.
## doi:10.18637/jss.v014.i06 <https://doi.org/10.18637/jss.v014.i06>.
##   - Zeileis A, Hothorn T (2002). "Diagnostic Checking in Regression
## Relationships." _R News_, *2*(3), 7-10.
## <https://CRAN.R-project.org/doc/Rnews/>.
##   - Zhu H (2021). _kableExtra: Construct Complex Table with 'kable' and Pipe
## Syntax_. R package version 1.3.4,
## <https://CRAN.R-project.org/package=kableExtra>.